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ABSTRACT 

A  discrete  potential  element  approach  to  subsonic  numer- 
ical lifting  surface  theory  has  been  developed  and  shown  to 
be  practical  in  predicting  the  nonsteady  loading  on 
harmonically  oscillating,  medium  to  low  aspect  ratio  wings. 
A  unique  method  of  including  the  wake  effect  in  the  wing 
kernel  function  matrix  prior  to  solution  of  the  singular 
integral  downwash  equation  was  devised,  thus  greatly 
simplifying  the  velocity  potential  formulation.   In  addi- 
tion, termination  of  the  effective  wake  a  finite  distance 
downstream  of  the  wing  was  investigated,  with  wing  loading 
found  to  converge  to  within  one  per  cent  in  an  effective 

4 

wake  length  of  four  root  chords. 

This  discrete  element  method  has  also  been  extended  to 
the  case  of  an  oscillating  wing,  cantilevered  from  a  cylin- 
drical fuselage,  to  investigate  nonplanar  interference 
effects.   This  interference  in  wing  loading,  while  of  rela- 
tively small  magnitude,  does  exist  in  both  pressure  amplitude 
and  phase  angle  distributions,  and  is,  therefore,  of  impor- 
tance in  three-dimensional  stability  analyses  of  wing/body 
configurations . 
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I.   INTRODUCTION 

The  aerodynamic  analysis  of  the  forces  over  an  oscillat- 
ing wing,  situated  in  a  steady  flow,  is  the  basic  problem 
considered  by  nonsteady  lifting  surface  theory.   The  deter- 
mination of  these  nonsteady  pressure  distributions  is  of 
prime  importance  in  the  consideration  of  stability.   A  basic 
objective  of  flutter  analysis,  for  example,  is  to  determine 
the  aerodynamic  loading  on  a  wing  oscillating  with  small 
amplitudes  in  a  definite  mode  and  with  a  given  frequency. 
Thus  the  problem  is  formulated  on  the  basis  that  the  motion 
and  deformation  of  the  lifting  surface  are  either  known  or 
composed  of  known  elemental  modes,  and  is  amenable  to  har- 
monic  modal  analysis.   A  downwash  integral  equation,  based 
on  the  transformed  wave  equation  for  linearized  potential 
flow,  relates  loading  (either  in  terms  of  velocity  potential 
or  pressure)  over  the  oscillating  wing  to  the  normal  velocity 
induced  by  the  wing  motion,  and  must  be  solved  numerically 
to  obtain  this  wing  loading. 

The  analysis  of  this  paper  will  deal  primarily  with  the 
subsonic  case,  this  being  the  one  in  which  the  largest 
number  of  applications  of  lifting  surface  theory  to  engineer- 
ing problems  are  to  be  found  [1] .   Most  of  the  work  in  sub- 
sonic nonsteady  lifting  surface  theory  is  based  on  the 
development  of  the  general  kernel  function  approach  by 
Kussner  in  19^0  [2],  which  was  first  formulated  successfully 
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for  the  computer  by  Watkins ,  Runyan,  and  Woolston  in  1955 
[31-   More  recent  techniques  have  been  developed  by  Stark  in 
1964  [4]  and  Laschka  in  1963  [5].   Solution  of  the  singular 
downwash  integral  equation  by  the  usual  lifting  surface 
methods  is  accomplished  through  assuming  the  loading  to  be 
a  series  of  preselected  functions  with  unknown  coefficients, 
and  then  determining  these  by  satisfying  the  normal  velocity 
condition  with  some  sort  of  collocation  technique.   The 
downwash  and  the  loading  are  related  by  kernel  functions, 
which  are  themselves  singular,  and  which  must  be  numerically 
evaluated. 

In  spite  of  the  fact  that  lifting  surface  theory  has 
enjoyed  a  great  deal  of  success,  certain  shortcomings, 
discussed  in  Refs.  1,  6,  and  7S  do  exist: 

(i)   Lifting  surface  theory  is  not  as  computationally 
economic  as  steady  or  two-dimensional  methods. 

(ii)   Wing  loading  functions,  which  incorporate  the 
proper  singular  behavior,  have  to  be  assumed  a  priori. 

(iii)   Uncertainty  exists  as  to  the  method  of  locating 
collocation  points. 

(iv)   Techniques  for  analyzing  control  surface  effects 
have  not  been  adequately  developed,  partly  because  of  (ii). 

(v)   Lifting  surface  theory  has  been  primarily  restricted 
to  planar  surfaces . 

In  order  to  attack  these  problems,  Houbolt  [6]  proposed  that 
the  wing  loading  be  represented  by  concentrated  pressure 
loads  in  the  form  of  Dirac  delta  functions.   Using  this 
loading  representation,  the  kernel  function  and  downwash 
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equation  were  formulated  and  analyzed  for  certain  restricted 
cases.   However,  no  general  numerical  solution  of  the  lifting 
surface   problem  was  developed. 

The  present  study  is  concerned  with  developing  such  a 
method  of  solution  to  the  subsonic,  nonsteady  lifting  surface 
problem,  in  terms  of  discrete  loading  elements.   However  in 
this  case  wing  loading  is  expressed  in  terms  of  velocity 
potential,  rather  than  pressure,  since  this  formulation 
allows  a  more  direct  approach  to  nonplanar  configurations. 
Haviland  [8]  reports  using  a  similar  approach  to  the  steady 
planar  case  with  consistent  results. 

The  aim  of  this  discrete  potential  element  method  is  to 
provide  a  new  approach  to  nonsteady,  subsonic  lifting 
surface  theory  which  eliminates,  or  minimizes,  the  problems 
previously  indicated.   Theoretical  development  of  this 
approach  was  made  from  basic  principles,  and  the  resulting 
formulation  applied  in  a  computer  program  which  analyzes 
harmonically  oscillating  wings  in  subsonic  flow.   A  unique 
way  of  including  the  wake  effect  was  also  developed  which 
greatly  simplified  the  velocity  potential  solution  of  the 
numerical  problem.   As  an  extension  of  the  theory,  and  an 
indication  of  the  versatility  of  this  discrete  element 
approach,  a  second  computer  program  was  developed  to  analyze 
interference  effects  on  oscillating  wings  mid-mounted  from 
steady  bodies.   To  the  author's  knowledge,  this  type  of 
wing/body  interference  investigation  has  not  previously  been 
presented  for  the  nonsteady  case. 


16 


II.   GENERAL  PROBLEM  DEVELOPMENT 

Following  the  method  of  development  by  Garrlck  [9],  the 
governing  equations  and  boundary  conditions  are  linearized 
by  application  of  small  perturbation  theory. 

A.   POTENTIAL  FLOW  EQUATIONS 

The  governing  equations  for  adiabatic,  irrotational, 
inviscid  flow,  as  developed  in  Ref.  10  are: 

Continuity 

2£  +  V  •  ptj  =  0  (2.1) 

at 

Momentum 

4 

dl  P 

Energy 

p 
—  =  constant  (y  the  specific  heat  ratio)       (2.3) 

PY 

Since  the  flow  is  irrotational,  the  curl  of  the  velocity 
vector  vanishes  (VxU=0)  and  a  velocity  potential  function 
can  be  defined,  such  that 

U  =  V$  (2.4) 

If  further,  a  perfect  gas  is  assumed,  then  from  the  energy 
equation  the  well  known  relationship 

8P/8p  =  a2  (2.5) 

is  obtained  for  constant  entropy  flows. 
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The  continuity  and  momentum  equations  can  now  be  rewritten  as 

1  Dp 

p  Dt 


^+V2$=0  (2.6) 


d  fa"*" 

where  =r  is  the  substantial  derivative  defined  by  N~r  +  U  •  V 

Combining  equations  (2.6)  and  (2.7),  a  single  equation 
in  $  is  obtained. 

or 

V2*-^^!*   0  (2.8) 

a  Dt 

This  is  the  governing  partial  differential  -equation  for 
general  nonsteady,  nonviscous,  potential  flow  of  a  perfect 
gas.   Equation  (2.8)  is  not  limited  to  small  disturbances, 
and  has  the  form  of  the  classic  wave  equation  in  terms  of 
the  substantial,  or  convective  derivative.   Thus,  the  flow 
disturbance  represented  by  the  velocity  potential  is  con- 
vected  by  the  local  fluid,  or  stream  velocity,  and  propogated 
as  a  wave  which  spreads  at  a  rate  equal  to  the  local  speed 
of  sound. 

In  this  general  form,  equation  (2.8)  is  highly  nonlinear 
due  both  to  the  quadratic  terms  in  the  convective  operator 
U*V,  and  to  the  interdependence  of  the  velocity  potential 
and  the  local  speed  of  sound.   Boundary  conditions,  in 
general,  will  depend  on  the  location  and  form  of  the  moving 
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body,  on  discontinuities  along  streamlines,  on  shock  waves, 
and  on  flow  conditions  at  infinity.   There  exist  no  general 
methods  for  obtaining  solutions  to  equation  (2.8),  and  it 
is  necessary  to  consider  small  perturbations,  to  linearize, 
to  reduce  the  number  of  equations,  or  to  try  schemes  of 
successive  approximations  to  attack  the  general  problem. 

B.   SMALL  PERTURBATION  THEORY 

Most  of  the  theoretical  developments  in  aerodynamics, 
including  nonsteady  aerodynamics,  are  based  on  the  concept 
of  small  perturbations.   There  are  two  aspects  to  the  prob- 
lem of  flow  past  a  body  creating  small  disturbances  in  an 
undisturbed  mainstream,  linearization  of  the  governing 
differential  equations  and  further  linearization  of  the 
boundary  conditions.  * 

Linearization  of  the  governing  equations  requires  the 
problem  to  be  posed  with  a  compressible  fluid  in  a  uniform 
stream  flowing  with  velocity  U   in  the  positive  x  direction. 
Small  disturbances  in  the  flow  are  defined  by 

U  =  (U^+u) i  +  v J  +  wk  =  U^i  +  u 

p  =  Poo+p        *  =  <(>«,+  $ 

_  -*■  — 

P  =  p   +  p        with  u  =  V(j) 

The  perturbation  velocities  are  considered  small  with  respect 
to  U  ,  a  ,  and  U  -a  ,  thus  ruling  out  transonic  flow.   The 
perturbation  density  and  pressure  are  also  restricted,  so 

that 
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p/pOT      and     p/poo«l 

Considering   4>    as    the  perturbation   velocity   potential, 
the    governing  equation    (2.8)    becomes 

2 


^2-r  1         9     ,    n        3 

V   *   "   —     It   +    U~   87 

oo 


$    =    0 


(2.9) 


where  higher  order  terms  in  the  U*V  operator  have  been 
neglected.   This  linearized  governing  equation  may  now  be 
applied  to  flow  problems  within  the  initial  small  perturba- 
tion assumptions.   Although  equation  (2.9)  was  developed  for 
a  uniform  flow  in  the  positive  x  direction  relative  to  a 
space-fixed  coordinate  system,  it  also  applies  to  the  case 
of  a  wing  moving  with  uniform  velocity  U^  in  the  negative  x 
direction  with  the  coordinate  system  fixed  to  the  wing. 

Within  the  concept  of  small  perturbations,  the  momentum 
equation  becomes 


8u   TT   9u      1  _- 

8t     »  3x      p 


(2.10) 


Substituting  the  perturbation  velocity  potential  (u  =  V<J)) 
into  this  equation,  a  linear  relationship  between  potential 
and  pressure  is  obtained. 


P  =  -P. 


Jt   +  u~  £? 


(2.11) 


Therefore  p  satisfies  the  same  differential  equation  (2.9) 
as  does  <}>. 

If  harmonic  motion  were  considered,  the  perturbation 
variables  would  have  the  form 
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$  =  <J>(x,y,z)e    ,  p  =  p(x,y,z)e 
u  =  u(x,y , z)e    ,  etc. 
and  the  governing  equation  would  become 

V2<f>  -  -4  (UM  ^  +  io))24>  =  0  (2.12) 

with  the  relationship  between  pressure  and  velocity  potential 

p  =  -pJuco  a!  +  ia))*  (2-13) 

In  the  above  equations  the  time  dependence  e    has  been 
factored  out.   This  relationship  between  p  and  $    can  be 
integrated  to  give 

-lOJX/U  „  •      r  /tt 

♦(x,y,z)    =   -        p   u  pU,y,z)e  dC      (2.14) 


CD       CO  —00 


where  far  ahead  of  the  disturbance,  the  perturbation  poten- 
tial is  assumed  to  be  zero. 

C.   LINEARIZED  BOUNDARY  CONDITIONS 

The  small  disturbance  assumption,  resulting  in  the 
linearized  governing  differential  equation  (2.9)  »    implies 
small  deviations  from  the  uniform  flow.   Therefore  it  is 
also  necessary  to  linearize  the  boundary  conditions  for  the 
physical  problems  considered,  such  as  properly  oriented  thin 
wings  and  slender  bodies.   These  boundary  conditions  consist 
of  the  mathematical  formulation  of  several  physical  and  phe- 
nomenological  statements  which  are  normally  grouped  as 
follows : 
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(i)   Surface  boundary  conditions.   No  flow  can  pass 
through  the  wing  or  body  surface,  that  is  to  say  the  total 
flow  is  tangential  to  the  body  surface. 

(ii)   Edge  conditions.   Sufficient  viscosity  is  present 
in  the  "nonviscous"  flow  to  determine  the  flow  pattern  near 
sharp  edges.   Thus  in  subsonic  flow  the  well  known  Kutta 
Condition  can  be  applied  to  the  wing  trailing  edge,  requiring 
that  the  surface  pressure  difference  remains  continuous 
near,  and  vanishes  at,  the  trailing  edge.   A  similar  condi- 
tion should  be  satisfied  at  the  side  edges.   Through  the 
method  of  matched  asymptotic  expansions,  Landahl  [11]  veri- 
fied that  at  trailing  and  side  edges  the  pressure  difference 
will  approach  zero  with  infinite  slope  due  to  weak  singu- 
larities in  the  first  derivative.   Conditions  at  the  leading 
edge  have  not  been  fully  explored,  but  for  small  disturbances 
to  rounded  leading  edges  it  is  sufficient  to  require  that 
the  total  integrated  force  be  finite  and  the  singularities 
in  the  pressure  distribution  be  of  the  proper  order  [9] • 
This  latter  requirement  is  fulfilled  with  the  pressure 
varying  as  the  inverse  square  root  of  the  distance  downstream 
of  the  leading  edge  [12]. 

(iii)   Wake  conditions.   The  free  vorticity  shed  from 
the  trailing  edge  in  the  wake  is  such  that  its  circulation 
together  with  the  bound  circulation  vanishes  in  accordance 
with  the  Helmholtz-Kelvin  Theorem.   It  is  assumed  that  the 
shed  wake  remains  where  it  is  formed,  floats  without  mutual 
interference  along  streamlines,  and  forms  a  continuous  sheet 
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of  discontinuity  coplanar  with  the  wing  in  the  direction  of 
flight.   Edge  effects  and  roll-up  of  the  sheet  at  infinity- 
are  ignored. 

(iv)   Conditions  at  infinity.   The  flow  is  considered 
uniform  at  infinity,  and  perturbation  waves  are  required  to 
propagate  away  from  the  sources  of  disturbances  and  to 
behave  properly  at  infinity. 

In  addition  to  the  above,  for  supersonic  flow  account 
must  be  taken  of  zones  of  influence  and  dependence,  as  in 
the  method  of  characteristics.   In  supersonic  flow  there  are 
also  other  problems  which  do  not  arise  in  the  subsonic 
analysis,  such  as  the  difference  between  subsonic  and  super- 
sonic edges. 

To  formulate  the  main  boundary  conditions,  a  thin  wing 
is  considered  lying  in  the  x-y  plane,  creating  small  dis- 
turbances in  a  uniform  stream  U^  flowing  in  the  positive  x 
direction.   The  linearized  tangential  flow  condition  on  the 

wing  surface  z   takes  the  form 
°  s 

w(x,y,zs,t)  =  {^+  U„^  (2.15) 

The  normal  fluid  velocity  w  can  be  expanded  in  a  power  series 
about  z  =  0,  so  that 

w(x,y,zs,t)  = w(x,y,0,t)  +  (•—    zg 

1  h   wl     -2  fo    -,£n 

+  27  7T   n  zs  +  (2'l6) 

I  9z  / z  =  0 
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In  addition,  problems  of  camber  and  thickness  can  be  separ- 
ated from  wing  motion  and  angle  of  attack  in  the  linearized 
formulation.   Therefore,  neglecting  higher  order  terms  in 
equation  (2.16),  the  tangential  flow  condition  becomes 


8z       8z 
w(x,y,0,t)  .  _£+  u.^ 


(2.17) 


where  z  (x,y,0,t)  represents  the  mean  camber  surface  posi- 
tion of  the  wing.  ' 

Since  w  is  an  even  function  of  z  (its  value  is  unchanged 
for  the  top  or  the  bottom  side  of  the  mean  surface),  the 
velocity  potential  <J>  is  an  odd  function  of  z,  as  is  the 
perturbation  pressure  p  through  the  linear  relationship  of 
equation  (2.11).   Therefore,  in  the  plane  of  the  wing  $ 
equals  zero  outside  of  the  wing  and  wake.   4>  does  not  equal 
zero  in  the  wake  because  of  the  discontinuity  in  the  u 
component  of  the  perturbation  velocity  across  the  wake.   The 
perturbation  pressure  is  also  zero  in  the  x-y  plane,  except 
at  the  lifting  surface  where  the  pressure  difference  Ap 
between  the  top  and  bottom  surfaces  is  given  by 


Ap(x,y,t)  =  p(x,y,0-,t)  -  p(x,y,0+,t) 
=  2p(x,y,0,t)  =  -2p 


$  (2.18) 


Figure  1  summarizes  the  formulation  of  these  boundary  condi- 
tions for  a  lifting  surface  in  the  x-y  plane. 

These  linearized  boundary  conditions  along  with  the 
linearized  governing  differential  equation  (2.9)  can  be  used 
to  attack  a  wide  range  of  perturbation  problems  in  subsonic 
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and  supersonic  flow.  Miles  [131  has  shown  that  a  sufficient 
condition  for  linearization  of  the  equations  for  the  pertur- 
bation potential  is  any  one  or  combination  of  the  factors 

|M2  -  1 I  »  0  (62/3) 

k  »  0  (62/3) 

M>>0(sV3) 

where  6,  M  6,  k6  ,  and  kM  6  must  be  sufficiently  small. 

3  00   '  '  00  " 

The   perturbation   equation   is   nonlinear   only   when   the 
following   conditions    are   jointly    satisfied 

|M2   -   1 I    =    0    (62/3) 

I         00  I  v  ' 

k   =    0    (62/3) 

AR   "    °    (6         } 

6  is  here  intended  to  be  a  thickness,  angle  of  attack,  camber, 
or  motion  parameter  serving  as  a  measure  of  the  flow 
disturbance . 
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III.   LIFTING  SURFACE  THEORY 

A.   BASIC  APPROACH 

The  basic  nonsteady  problem  considered  in  lifting  surface 
theory  is  that  of  a  linearized  thin  wing  harmonically 
oscillating  in  the  normal  direction  in  some  prescribed 
deflection  mode,  creating  small  disturbances  in  the  uniform 
mainstream.   Particular  solutions  to  the  linearized  governing 
differential  equation  (2.9)  are  required  which  give  pertur- 
bation velocity  potential  or  perturbation  pressure  distri- 
butions on  the  wing,  and  which  satisfy  the  essential  boundary 
conditions  discussed  in  Section  III-C.   Such  a  direct 
solution  has  been  possible  for  only  a  few  special  cases, 
such  as  two-dimensional  incompressible  flow  as  in  Refs.  14 
and  15.   For  the  more  general  three-dimensional  case,  an 
integral  equation  formulation  must  be  employed  which  makes 
use  of  a  downwash  kernel  function,  analogous  to  an  influence 
coefficient  development. 

Considering  the  time  dependence  e     in  the  harmonic 
analysis  as  being  factored  out,  the  integral  equation  relating 
the  downwash  amplitude  to  the  wing  loading  has  the  form 

w(x,y,0)  =  //  LU,n,0)K(x-£,  y-n,0)d?dn         (3-D 
S 

The  loading  L  can  either  be  the  velocity  potential  or  the 
pressure  distribution  over  the  surface,  and  K  is  the  kernel 

function  which  denotes  the  normal  velocity  (downwash)  on  the 
wing  at  (x,y,0)  due  to  a  unit  acoustic  singularity  located 
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at  (£,n,0).   In  other  words,  the  surface  is  modeled  by  a 
distribution  of  radiators,  formulated  in  terms  of  either 
potential  or  pressure  which  satisfy  the  linearized  governing 
equation 

V2<f> \   (U  ~  +  iuO2  (J)  =  0  (2.12) 

T    2   »  3x 
a 

oo 

w  is  the  downwash  at  point  (x,y,0)  on  the  wing,  determined 
by  the  harmonic  motion  of  the  surface  in  fulfilling  the 
tangential  flow  requirement.   All  of  these  quantities  are 
complex  in  the  spatial  domain  due  to  the  nonsteady  nature  of 
the  problem. 

Employing  the  velocity  potential  approach,  the  kernel 
function  is  defined  by 

K  =  w(x,y,0)  =  3<J>/3z  4  (3-2) 

where  <p  is  the  solution  to  equation  (2.12)  for  an  acoustic 
radiator  of  unit  strength.  In  this  formulation,  the  inte- 
gration (3-1)  must  be  carried  out  over  the  surface  of  both 
wing  and  wake,  since  the  potential  does  not  go  to  zero  in 
the  wake  region.  This  complicates  the  integration  process 
and  has  proven  a  major  disadvantage  of  the  potential  approach 

If  the  wing  is  modeled  by  a  pressure  distribution,  the 
kernel  function  is  denoted  by 

-iwx 


U-  x  U 


iu>£ 


K  =  w(x,y,0)  =  —     Lim        -^  f      p(£,y,z)e  d£ 

Poooo  z  +  0  -°°  (^^) 


from  equation  (2.14)  where  p  is  now  the  solution  to  equation 
(2.12)  for  an  acoustic  radiator  of  unit  strength.   The 
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integration  here  need  only  be  carried  out  over  the  wing 
surface,  since  the  pressure  goes  to  zero  in  the  wake. 
Because  of  this,  methods  based  on  the  pressure  formulation 
have  been  more  widely  employed.   This  kernel  function  (3-3) 
is  more  complicated  than  that  developed  from  equation  (3-2) 
and  has  not  been  integrated  in  closed  form  due  to  singulari- 
ties at  (y-ri)  =  0  and  (x-£)>_  0  .  However  these  singularities 
have  been  isolated  and  expressed  in  forms  which  can  be 
handled  by  numerical  procedures  [3] ■ 

The  acoustic  singularities  used  to  model  the  wing  surface 
in  subsonic  flow  are  harmonically  pulsating  doublets  (dipoles) 
oriented  in  the  z,  or  wing-normal,  direction.   Doublets  are 
employed  in  order  to  represent  the  pressure  difference 
between  the  two  surfaces  of  the  wing,  while  in  supersonic 
flow  acoustic  sources,  or  monopoles,  are  used  due  to  the 
independence  of  the  two  wing  surfaces  in  that  flow  regime. 
Kiissner  [2]  developed  the  general  formulation  of  the  pressure 
doublet  kernel  function,  through  use  of  the  Lorentz  trans- 
formation, in  order  to  apply  solutions  to  the  classical  wave 
equation  of  physics 

V2$   1  ■  A  .  0 

a        9t 

to  the  governing  differential  equation  of  lifting  surface 
theory  (2.9).   This  method  of  attack  is  used  in  Appendix  A 
to  develop  the  kernel  functions  used  in  this  report. 
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B.   NUMERICAL  SOLUTIONS  TO  LIFTING  SURFACE  THEORY 

To  find  the  potential  or  pressure  loading  on  the  lifting 
surface  actually  requires  the  inverse  solution  of  equation 
(3.1) »  since  the  downwash  velocity  distribution  is  deter- 
mined by  the  type  of  surface  motion  prescribed  through 
equation  (2.15).   The  numerical  methods  that  have  been 
developed  in  subsonic  lifting  surface  theory  rest  on  this 
inverse  solution  to  the  singular  integral  equation.   The 
usual  procedure  is  to  assume  the  loading  to  be  a  series  of 
preselected  functions  with  unknown  coefficients  and  then 
to  determine  these  by  satisfying  the  downwash  velocity 
distribution  exactly  in  a  set  of  collocation  points  [16,  17], 
or  approximately  in  the  least  squares  sense  in  a  larger  set 
of  control  points  [43  18]  ,  or  by  satisfying  certain  integral 
relations  derived  from  variational  procedures  [19,  20]. 

Ashley  and  Landahl  [21] ,  and  Landahl  and  Stark  [1]  have 
presented  survey  papers  summarizing  the  status  of  numerical 
lifting  surface  theory,  and  presenting  formulations  of  the 
integral  downwash  equation  and  the  kernel  function  in  the 
various  methods  of  attack,  such  as: 

(i)   Velocity  potential.   As  previously  discussed,  the 
kernel  function  used  in  this  formulation,  while  singular,  is 
much  simpler  than  that  employed  with  the  pressure  loading 
approach.   However,  integration  must  be  carried  out  over 
both  wing  and  wake.   The  velocity  potential  formulation  has 
been  employed  by  Jones  [22]  for  M^  =  0  in  1952. 

(ii)  Acceleration  potential.   Formulating  this  approach 
in  terms  of  an  acceleration  potential  defined  by  c  =  —  2 \p  [23], 
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allows  the  use  of  the  same  kernel  function  as  in  the  velocity 
potential  formulation,  but  removes  the  necessity  of  integrat- 
ing over  the  wake,  since  c   =0  there.   However,  the  solution 
is  not  unique,  since  multiples  of  eigen  solutions  with 
J  =  0  on  the  wing  may  be  added,  and  integration  must  be 

dZ 

extended  into  the  region  ahead  of  the  leading  edge  to  achieve 
uniqueness.   This  formulation  has  been  used  for  two- 
dimensional  cases  and  also  for  certain  three-dimensional 
cases  [24]  . 

(iii)   Pressure  formulation.   The  original  development 
by  Kussner  [2]  has  been  formulated  for  the  computer  by 
Watkins,  et  al  [3]  and  by  Richardson  [16].   This  procedure 
provides  direct  determination  of  the  pressure  with  integra- 
tion required  only  over  the  wing  surface,  but  the  kernel 

4 

function  is  highly  singular  and  needs  to  be  evaluated  through 
numerical  quadrature,  thereby  increasing  computer  time 
considerably . 

Various  refinements  of  the  basic  approaches  discussed 
above  have  been  developed,  such  as  the  integrated  acceler- 
ation potential  by  Stark  [4]  and  the  advanced  velocity 
potential  [1] ,  but  these  do  not  change  the  basic  problem 
formulation,  nor  the  basic  advantages  and  disadvantages  of 
the  continuous  loading  function  methods. 

In  lifting  surface  theory  the  loading  is  expressed  as 
a  linear  combination  of  preselected  functions,  with  coeffi- 
cients to  be  determined  by  satisfying  the  tangential  flow 
condition  through  the  integral  equation  (4.1).   In  selecting 
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the  loading  functions,  known  results  from  two-dimensional 
incompressible  flow  theory  have  been  used.   The  loading  is 
separated  into  chordwise  and  spanwise  functions,  so  that  a 
pressure  function  element  would  be  represented  as 

p   =  f  (C)g  (n) 
^mn    m   ton 

These  loading  functions  have  been  represented  by'  Fourier 
expansions  [17],  power  series  [25,  26],  Tschebycheff 
polynomials  [27,  28],  Legendre  polynomials  [4],  and  so  on. 
The  basic  requirement  is  that  these  functions  themselves 
satisfy  the  applicable  boundary  conditions  from  Section  II, 
such  as  leading,  side,  and  trailing  edge  behavior. 

C.   SOLUTIONS  BASED  ON  DISCRETE  LOADING  LINES 

In  steady  flow,  the  lifting  line  theory*  of  Prandtl 
employed  a  discrete  vortex  lifting  line,  with  variable 
strength,  placed  at  the  quarter  chord  line.   Wieghardt  [29] 
extended  this  method  for  rectangular  wings  by  using  several 
lifting  lines,  while  Rubbert  [30],  Dulmovits  [31],  and 
Hedman  [32]  assumed  constant  strengths  of  the  lifting  lines 
in  subintervals  in  a  vortex  lattice  approach. 

The  vortex  lifting  line  in  steady  flow  corresponds  in 
the  nonsteady  case  to  a  potential  doublet  strip  with  strength 
varying  harmonically  in  the  streamwise  direction.   The  use 
of  a  discrete  lifting  line  approach  to  the  nonsteady  case 
was  first  proposed  by  Jones  [22]  for  the  incompressible 
case,  and  by  Runyan  and  V/oolston  [25]  for  general  subsonic 
flow. 
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As  in  the  steady  case,  this  lifting  line  approach  for 
nonsteady  motion  has  been  extended  to  the  doublet  lattice 
method  by  Stark  [1] ,  and  by  Albano  and  Rodden  [331.   The 
downwash  from  the  doublet  lattice  strip  is  closely  related 
to  the  kernel  function  discussed  previously  in  this  section. 
Since  the  surface  loading  is  replaced  by  discrete  doublet 
strips,  it  is  not  necessary  to  assume  the  loading  functions, 
as  required  in  the  continuous  function  approach  to  lifting 
surface  theory,  and  the  integral  downwash  equation  (3.1)  is 
replaced  by  a  matrix  equation  relating  downwash  on  the  wing 
to  doublet  lattice  strength.   This  approach  then  removes 
one  of  the  shortcomings  to  general  lifting  surface  theory, 
but  still  retains  the  remaining  features. 

D.   NONPLANAR  CONFIGURATIONS 

Limited  work  has  been  done  in  the  area  of  general  non- 
planar  configurations,  since  interference  investigations 
have  been  primarily  involved  with  combinations  of  lifting 
planar  surfaces.   The  kernel  functions  for  interfering  planar 
surfaces,  as  developed  by  Davies  [3^]  and  by  Landahl  [35]  s 
become  more  complicated  but  do  not  contain  any  new  singulari- 
ties.  In  the  selection  of  suitable  loading  functions,  the 
behavior  at  the  intersection  of  lifting  surfaces  must  be 
taken  into  account.   The  T-tail  has  been  treated  by  Stark 
[^]  and  by  Davies  [35]  }  while  calculations  for  a  delta  wing 
with  folded  tips  have  been  given  by  Vivian  and  Andrew  [36]. 
Lashka  [5]  has  analyzed  the  effects  of  wing  tip  pylons  as 
well  as  interfering  planar  surfaces  [37],  while  Rodden  [38] 
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has  included  wing/body  and  wing/pylon  effects  in  nonsteady 
wing  loading  in  subsonic  flow. 
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IV.   DISCRETE  POTENTIAL  ELEMENT  DEVELOPMENT 

A.   WING  ANALYSIS 

1 .   Basic  Formulation 

The  integral  equation  relating  downwash  to  wing 
loading  (either  velocity  potential  or  pressure)  in  lifting 
surface  theory,  as  discussed  in  Section  III,  is 

w(x,y,0)  =  //  L(£,n,0)K(x-c,y-n,0)d€dn        (3-D 
S 

where  harmonic  motion  is  assumed  so  that 

w(x,y,z,t)  =  w(x,y,z)ela) 
L(x,y,z,t)  =  L(x,y,z)elw 

4 

In  the  subsequent  development,  the  time  dependence  will  be 
considered  factored  out,  so  that  all  symbols  will  refer  to 
complex  valued  spatial  variables.   The  downwash  w  is  deter- 
mined by  the  boundary  condition  of  no  flow  through  the  wing 
surface  as  given  for  general  motion  by  equation  (2.7),  and 
here  for  harmonic  motion 

dz 

w  =  U   t-^-  +  iwz  (4.1) 

00  8x       s 

where  z   is  the  amplitude  of  the  surface  motion  prescribed 
for  the  problem.   The  form  of  acoustic  radiator  used  to 
model  the  wing  is  taken  as  a  potential  doublet,  with  axis  in 
the  z,  or  wing  normal,  direction.   This  dipole  singularity 
is  employed  because  of  the  requirement  to  sustain  a  pressure 
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difference  across  the  wing  in  subsonic  flow.   K  is  the  kernel 
function  relating  the  downwash  at  point  (x,y,0)  on  the  wing 
due  to  a  unit  potential  doublet  located  at  (£,n,0). 

In  normal  lifting  surface  theory,  the  loading  L  would  be 
represented  by  continuous  chordwise  and  spanwise  functions 
which  satisfy  the  surface  boundary  conditions  developed  in 
Section  II.   However,  in  the  discrete  potential  element 
approach  this  loading  is  represented  by  a  network  of  Dirac 
delta  functions  formulated  in  terms  of  the  perturbation 
velocity  potential.   Thus  the  velocity  potential  distribution 
over  the  wing  is  replaced  by  a  series  of  point  functions  in 
the  form  of  harmonically  oscillating  potential  doublets  with 
axes  in  the  z  direction.   In  addition,  the  wake  must  also 
contain  a  network  of  these  doublets  since  the  velocity 

4 

potential  Is  not  zero  there  (Section  II)  and  the  integration 
region  of  equation  (3-1)  must  include  both  wing  and  wake.   A 
representation  of  this  potential  doublet  grid  is  shown  in 
Figure  2. 

In  the  discrete  formulation,  the  integral  equation  (3-1) 
is  replaced  by  the  matrix  equation 

{w}  =  [K]  {4,}  (4.2) 

where  w.  is  the  downwash  at  the  i    point  on  the  wing,  and 

$.  is  the  strength  of  the  j    potential  doublet  on  the  wing 
j 

and  wake  in  the  grid  of  Figure  2.   The  location  of  the  down- 
wash  control  points  are  coincident  with  the  potential  doublet 
positions  at  the  center  of  each  wing  control  box.   The  choice 
of  this  location  will  be  discussed  in  Section  IV. A. 4.   In- 
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modeling  the  wing,  the  velocity  potential  and  the  pressure 
are  considered  constant  over  each  of  the  control  boxes,  with 
approximations  to  continuous  distributions  given  by  these 
values  at  the  control  box  midpoints. 

The  kernel  function  matrix  elements  represent  the  down- 
wash  due  to  potential  doublets  of  unit  strength  which  satisfy 
the  governing  linearized  differential  equation 

V24>    -  -\    (UM    3^  +    ico)24>    =    0  (2.12) 

a 

00 

As  developed  in  Appendix  A,  this  kernel  function  has  the 
form 

K,,  =  ^-~  (1+i  — ^R)exp{i  — SL.  [M(x,-x.)-R]}(4.3) 
1J    4ttRJ      a  ^  a  r      X     3 


00  00 


where 


/(x.-x  )2  +  32(y1-y1) 


In  this  formulation  the  kernel  function  is  singular  only 
when  the  downwash  control  point  is  coincident  with  the  poten- 
tial doublet  location  (i=j).   The  handling  of  this  singular- 
ity is  vital  to  this  development  and  is  discussed  in  detail 
in  Section  IV. A. 4. 

The  discrete  loading  element  approach  to  the  nonsteady 
wing  problem  requires  the  solution  of  the  system  of  linear 
complex  equations  (4.2)  to  determine  the  velocity  potential 
vector  {<{)}.   The  downwash  vector  {w}  is  specified  by  the 
wing  motion  through  equation  (4.1),  so  that  the  solution  has 
the  form 
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{<}>}    =     [K]"1    {w}  (l|.i|) 

where  {<t>}  and  {w}  cover  both  wing  and  wake  regions.   The 
resultant  velocity  potential  distribution  is  used  to  deter- 
mine the  pressure  loading  over  the  wing  through  the  rela- 
tionship 

Ap  =  -2poo(Uoo  g|  +  io>)4>  (4.5) 

taken  for  the  harmonic  case  from  equation  (2.18). 

2.   Wake  Effect 

One  of  the  historic  drawbacks  to  the  potential 
approach  to  lifting  surface  theory  has  been  the  necessity 
of  including  the  wake  region  in  the  solution  of  the  integral 
equation  (3.1).   This  equally  complicates  the  discrete 
element  approach  since  the  matrix  equation *( 4 . 2 )  is  required 
to  include  the  effect  of  the  wake  potential  dipole  grid. 
Even  though  the  effective  wake  is  terminated  at  some  repre- 
sentative length  behind  the  wing,  as  is  commonly  done  in 
steady  flow  problems,  the  size  of  the  resulting  potential 
strength  vector  and  kernel  function  matrix  would  severely 
restrict,  or  entirely  preclude,  the  use  of  this  method  even 
on  large  modern  computers. 

Following  a  suggestion  by  Professor  R.  E.  Ball  of  the 
Naval  Postgraduate  School,   a  careful  examination  of  the 
boundary  conditions  governing  harmonic  wing  motion,  as 
depicted  in  Figure  3>  was  made.   In  the  wake  region  the 
pressure  is  zero,  but  the  velocity  potential  is  not  zero 
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because  of  the  discontinuity  in  the  u  component  of  the  per- 
turbation velocity.   <J>  and  p  are  related  by 

In  the  wake,  therefore,  this  relationship  becomes 

Uw  §±  +  lu+  =  0  (4.5) 

or 

34      .   a)  . 

87  ~  "  1  U"  * 

00 

which  can  be  integrated  to  give 

-ifT  (xw-xm} 
*w  =  V  e  (4.6) 

<}>   is  the  strength  of  the  wake  velocity  potential  at 

(x  ,y.,0),  and  <j>   is  the  strength  of  the  ve'locity  potential 
w   i  m 

at  the  boundary  between  the  wing  and  the  wake  (x  ,y.,0) 
along  the  line  of  integration  y=y..  In  the  discrete  poten- 
tial element  formulation,  d>   becomes  the  strength  of  the 

'   m  ° 

last  wing  dipole  on  the  appropriate  chord  line,  with  (x  -x  ) 

w   m 

the  distance  between  this  dipole  and  the  wake  dipole  <J>  , 

w 

as  is  shown  schematically  in  Figure  4. 

The  wake  potential  strength  distribution  is  therefore 
expressible  as  a  function  of  the  wing  potential  distribution 
and  can  be  included  in  the  wing  doublet  kernel  functions. 
This  relationship  can  be  represented  by 

Uw>  =  [A]  {<*>} 
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where  {<J>}  is  the  wing  velocity  potential  vector  and  the 
appropriate  elements  of  [A]  express  the  relationship  of 
equation  (4.6)  with  respect  to  the  wake  potential  vector 
{<p    }.   The  downwash  matrix  equation  (4.2)  therefore  becomes 

{w}  =  [Kx]  {<*>}  +  [K2]  Uw) 

=  [^1  W    +    [K2]  [A]  {<(>} 

=  [K1+K2A]  {4.} 

The  size  of  this  matrix  equation  has  thus  been  limited  to 
that  necessary  to  represent  only  the  wing  grid,  making  the 
discrete  potential  element  approach  much  more  tractable  on 
the  computer. 

3.   Matrix  Equation  for  Symmetric  Motion 

In  considering  a  wing  undergoing  harmonic  oscilla- 
tions which  are  symmetric  with  respect  to  the  x-z  plane,  the 
wing/wake  planform  can  be  represented  as  shown  in  Figure  5« 
The  x-z  plane,  being  a  plane  of  symmetry,  represents  a  no- 
flow  surface,  or  reflection  plane,  for  the  disturbances 
caused  by  the  symmetric  motion  of  the  full-span  wing.   View- 
ing the  physical  problem  in  a  different,  but  equally  valid 
light,  the  half-span  wingQ,  with  its  root  along  the  x  axis, 
can  be  considered  undergoing  oscillations  in  the  presence  of 
an  infinite  wall  coincident  with  the  x-z  plane.   Therefore, 
the  motion  of  the  wing  need  only  be  prescribed  for  wing  Q, 
while  the  function  of  the  virtual  wing  Q  and  wake  ©  is  to 
establish  the  no-flow  'condition  at  the  x-z  plane. 

The  full  wing/wake  system  is  represented  by  the  matrix 
equation 
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where  subscripts  denote  the  areas  of  Figure  5  and  K. .  is  the 
kernel  function  matrix  for  downwash  at  points  in  area  ©  due 
to  doublets  located  in  areaQ.   For  symmetric  motion 

U  }  =  (4^}    and   {c^}  =  {<J>2> 
so  the  required  downwash  equation  is  reduced  to 


{wn}  =  [Knn+Kno   Kno+K-,„]  <--- 


11  "13     12   14 


*2 


In  the  previous  section,  the  relationship  between  wing  and 
wake  potential  distributions  was  developed,  and  can  here  be 
expressed  as 

{4>2>  =  [A]  {<|)1} 

Therefore,  the  final  wing  downwash  matrix  equation  has  the 
form  of  equation  (4.2) 


{w1>  =  [K]  {<j>1} 


but  the  kernel  function  matrix  is  here  formed  from  the 
following  expression 

[K]  =  [Kni]  +  [K,  ,]  +  [Kno+KnJl]  [A] 


(4.7) 


'11 


'13 


12TIV14 


(4.8) 


Thus,  through  use  of  motion  symmetry  and  the  boundary  condi- 
tion in  the  wake,  the  effective  integration  of  equation  (4.7) 
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need  only  be  carried  out  over  the  half-span  wing  surface  Q. 
Acknowledgement  of  the  full  physical  problem  is  obtained 
through  formation  of  the  kernel  function  matrix  via  equation 
(4.8),  which  incorporates  effects  of  the  potential  dipole 
grids  in  all  four  areas  of  the  wing  and  wake  depicted  in 
Figure  5. 

It  should  be  noted  that  antisymmetric  motion  of  the  wing 
about  the  x-z  plane  requires  that 

{$  }  =  -{01)    and   {<J>4>  =  -{<j>2> 

The  final  kernel  function  matrix  is  therefore  formed  by 
[K]  =  [K1:L]  -  [K13]  +  [K12-Kl4]  [A] 

Since  any  general  harmonic  motion  of  the  wing  can  be 
expressed  as  a  combination  of  symmetric  and*  antisymmetric 
modes,  this  approach  has  general  application  and  need  not 
be  restricted  to  the  symmetric  case  considered  in  detail 
here . 

4.   Doublet  Singularity 

The  location  of  the  collocation  points  required  to 
satisfy  the  no-flow  condition  in  lifting  surface  theory  is 
historically  based  on  two-dimensional  steady  flow  theory, 
Vortex  lattice  methods  such  as  Rodden's  [38]  arrange  the 
vortex  strip  on  the  local  quarter  chord  of  the  control  box, 
with  the  downwash  collocation  point  centered  on  the  local 
three-quarter  chord  line  as  in  two-dimensional  thin  airfoil 
theory.   Houbolt  [6]  proposes  using  this  local  quarter  chord, 
three-quarter  chord  control  box  grid  in  conjunction  with 
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concentrated  pressure  loads.   This  type  of  grid  network  was 
tried  by  the  author  in  an  early  form  of  the  wing  analysis 
computer  program  with  unsatisfactory  results.   The  approach 
was  found  not  to  converge,  but  to  be  very  sensitive  to  grid 
size;  that  is,  to  the  distance  between  the  potential  dipole 
and  its  associated  control  point.   This  is  much  like  the 
sensitivity  that  a  continuous  loading  method  experiences 
when  a  control  point  is  located  too  close  to  a  wing  edge 
[11]. 

The  control  or  collocation  points  were  subsequently 
placed  at  the  center  of  the  wing  control  boxes  (Figure  2) 
coincident  with  the  potential  doublet  locations.   In  this 
way,  the  above  mentioned  sensitivity  to  grid  spacing  was 
removed,  but  the  value  of  the  upwash  of  a  doublet  at  its 
own  control  point  had  to  be  determined.   As  can  be  seen 
from  equation  (4.3),  the  kernel  function  (K..)  is  singular 
at  the  doublet  location  (i=j).   This  singularity  in  the 
upwash  from  a  doublet  is  pictured  in  Figure  6a,  where  the 
doublet  produces  an  infinite  upwash  at  its  location,  but 
finite  and  decreasing  downwash  in  the  plane  of  the  wing  as 
the  distance  from  the  doublet  is  increased. 

If  the  dipole  is  considered  within  the  framework  of  the 
discrete  element  grid  where  dipole  strengths  and  downwash 
velocities  are  held  constant,  or  averaged,  over  each  control 
box,  a  cross-section  of  the  dipole  flow  pattern  would  appear 
as  in  Figure  6b,  where  the  infinite  upwash  at  the  doublet 
has  been  replaced  by  a  finite  value  w  .   To  determine  a 
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value  for  w  which  adequately  represents  the  singularity 
strength  within  the  framework  of  the  discrete  element 
approximation,  the  law  of  continuity  was  employed.   The 

total  fluid  outflow  in  the  discrete  representation  is  w  A  , 

o  o ' 

where  A   is  the  area  of  the  doublet  control  box.   Modeling 
the  entire  x-y  plane  with  control  boxes  without  regard  to 
wing/wake  geometry,  the  total  amount  of  fluid  passing  back 

00 

through  the  x-y  plane  is   I   w  A  ,  where  w   is  the  average 

n=1      n   n  n 

velocity   over   the    control   box  with   area  A    .      Therefore, 
j  n  ' 

continuity    requires 


wA+EwA=0 
o   o  ,      n   n 

n=l 

and  if  the  A  '  s  are  chosen  so  that  A   =  A  ,  the  upwash 
n  o    n'      ^ 

velocity  is  determined  by 


w  =  -   E   w  (4.9) 

o       ,   n 
n=l 

In  actual  practice,  the  summation  of  the  downwash  velocities 

through  the  x-y  plane  is  necessary  only  to  some  finite  radius 

from  the  doublet,  since  the  inverse  proportionality  of  the 

kernel  function  with  distance  from  the  doublet  causes  the 

value  of  w   to  converge  within  a  reasonable  summation. 
o  ° 

5 .   Section  and  Wing  Coefficients 

To  obtain  section  lift  and  moment  coefficients,  the 
chordwise  pressure  distribution  at  each  spanwise  station 
was  integrated  in  the  manner  employed  in  thin  airfoil  theory. 
First  a  coordinate  transformation  of  the  chordwise  variable 
is  made  as  follows 
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x/c  =  -^(i-cos  e) 


The  pressure  coefficient  is  then  expressed  in  terms  of  a 
Fourier  expansion  of  the  new  variable 

00 

c  (6)  =  c    cot  |  +   2   c   sin  n9  (4.9) 

p       po         n=l   pn 

where  the  first  term  accounts  for  the  leading  edge  singu- 
larity, while  the  trailing  edge  slope  singularity  is 
acknowledged  by  the  infinite  series.  The  section  lift  is 
obtained  from 


ct  -  /   cp  d[f|  (4.10) 

Substituting  the   Fourier  expansion  for   c      and  performing   the 

coordinate    transformation,    equation    (4.10)    becomes 
c  c 

P     TT  °°     p     77 

cp  =  — -  /   (l+cos0)d9  +   E   — -  /  sin  n9  sinGde 

2   0  n=l   2   0  (JL11) 

Due  to  the  orthogonality  of  the  sine  function,  this  integra- 
tion yields 

c0  =  \    (c    +  \   c   )  (4.12) 

*■    2    pQ    2   p1 

The  section  moment  about  the  leading  edge  is  defined  by 
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Performing  the  same  transformation  and  integration  as  for 
the  section  lift  coefficient,  the  following  results 

cmo  =  "  I   (cpo+  cPl"  2  cp2)  (4.14) 

Transferring  this  moment  coefficient  to  the  mid-chord  point 
requires  the  further  calculation 

c   =  c  +  ^  c0  (4.15) 

m-,     m   2  I 
h  o 

The  integration  of  section  lift  and  moment  values  to  wing 
coefficients  is  performed  in  a  similar  manner.   Here  the 
coordinate  transformation  is  made  to  the  spanwise  variable, 
such  that 

y/b  =  |  (1-cos  6) 

The  product  of  the  section  chord  and  the  section  lift 
coefficient  is  then  expressed  as  the  Fourier  series. 

oo 

a 

cc.,  =  c,  cos  —   +   I   c,  sin  n9  (4.16) 

o        n=l   n 

Performing  the  integration 

of  the  Fourier  expansion  of  the  section  lift  in  the  trans- 
formed coordinate  system,  the  following  results 

CT  =  £  U  c   +  ?  c   1  (4.18) 

L   s  I  3   A   T  A,  I 

\     o      1  / 
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The  same  relationship  is  obtained  for  the  wing  moment 
coefficient  when  the  sectional  moment  is  expressed  as  in 
the  Fourier  expansion  of  equation  (4.16).   For  swept  wings, 
the  sectional  moments  were  transfered  to  an  axis  extending 
from  the  half-root-chord  point,  prior  to  developing  the 
Fourier  series. 

The  infinite  series  of  equations  (4.9)  and  (4.16)  are 
of  course  terminated  at  the  number  of  chordwise  and  span- 
wise  points  respectively  of  the  wing  discrete  element  grid. 
This  places  a  minimum  limit  on  the  number  of  chordwise 
and  spanwise  control  points  which  are  required  to  achieve 
valid  sectional  and  wing  coefficient  values,  a  factor  which 
will  be  discussed  in  Section  VI. 

B.   WING/BODY  ANALYSIS 

1.   Basic  Formulation 

Inclusion  of  a  finite  radius  body  in  the  nonsteady 
lifting  surface  analysis  is  equivalent  to  adding  one  more 
boundary  condition  to  the  problem  formulation:   the  require- 
ment for  no  flow  through  the  body  surface.   For  the  analysis 
pursued  here,  the  body  will  be  stationary  with  regard  to 
the  perturbation  motion,  and  will  be  idealized  as  an 
infinitely  long  cylindrical  surface  with  axis  coincident 
with  the  undisturbed  flow.   The  cantilevered,  midmounted 
wing  is  harmonically  oscillating  within  the  limitations  of 
small  disturbance  theory  previously  developed. 

In  steady  flow  analysis,  the  body  effects  have  been 
traditionally  handled  by  singularities,  matching  the  wing 
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singularity  distribution,  at  image  points  within  the  body  in 
the  wing  plane  [391-   This  method  will  satisfy  the  nonsteady 
boundary  conditions  only  in  a  quasi-steady  sense,  because 
of  phase  differences  between  disturbances  at  the  body  sur- 
face caused  by  a  wing  singularity  and  its  corresponding 
image.   Therefore,  to  model  the  body  in  the  nonsteady 
problem,  a  system  of  singularities  are  placed  on  the  body 
surface  establishing  a  grid  similar  to  that  for  the  wing. 
These  curvilinear  panels  each  have  an  harmonically  oscillat- 
ing singularity  at  its  center,  coincident  with  a  normalwash 
control  point.   Thus,  in  effect  the  wing  grid  is  merely 
extended  over  the  body  surface  as  shown  in  Figure  7- 

In  order  to  handle  the  more  complex  geometry  of  the 
wing/body  problem,  a  coordinate  transformation  from  the 

4 

rectangular  (x,y,z)  system  of  the  wing  analysis  is  made  to 
a  cylindrical  (r,6,x)  system,  where  the  undisturbed  flow 
direction  (x  axis)  is  common.   Thus  the  transformation  is 
defined  by 

x  =  x 

y  =  r  cos  9  (4. 19) 

z  =  r  sin  6 

Figure  8  shows  the  wing/body  configuration  with  the  boundary 
conditions  in  the  cylindrical  coordinate  system.   The  down- 
wash  velocity  over  the  wing  is  now  related  to  the  velocity 
potential  distribution  by 
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and  the  no-flow  condition  for  the  body  surface  is  stated  as 


u 


8r 


(4.21) 


r=r 


B 


where  rR  is  the  body  radius.   As  in  the  wing  analysis,  the 
wing  downwash  velocity  ufi  is  determined  from  the  type  of 
surface  motion  prescribed  for  the  problem. 
2 .   Governing  Matrix  Equation 

The  governing  matrix  equation  for  the  wing/body 
problem  takes  on  a  much  more  formidable  appearance  than  the 
basic  equation  (4.2)  for  the  wing  along 
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(4.22) 


where    ^uo^  -   wing  downwash  vector 

{u  }  =  body  normalwash  vector 

{4>y}  =  wing  singularity  potential  vector 

{<j>g)  =  body  singularity  potential  vector 

The  kernel  function  matrices  are  defined  as  follows 

Ky  =    downwash  on  wing  due  to  wing  singularities 
KWR  =  downwash  on  wing  due  to  body  singularities 
KRW  =  normalwash  on  body  due  to  wing  singularities 
Kg  =   normalwash  on  body  due  to  body  singularities 

Direct  solution  of  equation  (4.22)  for  the  velocity  potential 
distribution  would  be.  limited  due  to  the  computer  storage 
requirements  of  the  large  complex  kernel  function  matrix  . 
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However  this  approach  is  not  feasible  at  any  rate,  since 
the  kernel  function  matrix  is  ill  conditioned  in  this  form 
and  does  not  lend  itself  to  efficient  inversion. 

To  attack  the  problem,  equation  (4.22)  is  separated  into 
two  matrix  equations 

<V  =  tKW]{*W}  +  1KWB1{*B} 
{ur}  =  IKBWH*W}  +  [KB]{*B) 

The  second  of  these  can  he  solved  for  the  body  potential 
distribution 

UB)  =  [Kg]"1  [ur-KBW*w]  (4.23) 

which  is  then  substituted  into  the  first  equation. 

<V  ■  [KWH*W}  +  [IW  [KBrl{V  -  [IW  [KBrl[KBW!  CV 
This  may  be  rewritten  in  the  following  form 

N-^raV1  ur}  =  'VWb1  kbwhV  <*-2*> 

The  left  hand  matrix  is  a  modified  wing  downwash  vector 
incorporating  body  effects,  as  does  the  modified  kernel 
function  matrix  on  the  right  hand  side  of  equation  (4.24). 
In  the  analysis  considered  here  the  body  is  steady  at 
zero  angle  of  attack  to  the  main  stream,  while  the  wing  is 
undergoing  nonsteady  motion,  so  that 

{ur>  =  {0} 

Therefore,  equation  (4.24)  becomes 
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{ue}  ■  [Wb1  kbw]{V  (4-25) 

which  is  in  the  same  form  as  equation  (4.2) a  but  in  which 
the  modified  kernel  function  matrix  incorporates  the  body 
boundary  condition. 
3.   Symmetry 

The  individual  kernel  function  matrices  of  equation 
(4.25)  incorporating  body  effects  can  be  quite  large  when 
control  points  are  placed  around  the  entire  circumference 
of  the  body,  severely  limiting  the  utility  of  this  method. 
However  these  matrices  may  be  reduced  appreciably  in  size 
through  use  of  symmetry.   Considering  the  representation  of 
Figure  9a,  it  can  be  seen  that  for  wing  motion  symmetric 
about  the  x-z  plane 

{*W2»  ■  {*Wl};   {*B2}  "  {+Bx};   {*B3}  "  {V  "-26) 

where  { $„   }  stands  for  the  doublet  strength  distribution  on 

i 
the  respective  wing,  and  {<J>R  }  stands  for  the  singularity 

i 
strength  distribution  on  the  body  in  the  respective  quadrant 

It  can  also  be  seen  that  due  to  the  antisymmetric  nature  of 

the  wing  doublet  flow  with  respect  to  the  x-y  plane 

U   }  =  -  {<f>   };   U   }  =  -  {(f)   }  (4.27) 

4         1       3         2 

since  the  purpose  of  the  body  singularities  is  to  counter 
the  flow  of  the  wing  doublets  through  the  body  surface. 

Considering  first  the  body  kernel  function  matrix  (Kg) , 
which  is  defined  by  the  relation 
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{ur>    =    [KB]UB} 


or 


{u    }    =    [K      ]  U      }  +   [KR    ]  UR    }  +  [KR    ]  {$      }  +   [Kn    ]  {(j)p    } 

r       B1    B1      B2    B2      B^    B^      B^    B^ 

Symmetry  allows  this  formulation  to  be  reduced  to 

{ur>  =  [KB  +KB  -Kfi  -KB  ]{*B  >  (4.28) 

In  the  case  of  the  wing-body  kernel  function  matrix  (KWR) , 
this  same  consideration  produces 

{V  ■  [KWB1+KWB2-KWB3-KWB4H*Bl}  °'29) 

Control  points  need  only  be  placed,  therefore,  over  the 
first  quadrant  of  the  body  surface,  reducing  the  size  of 
the  body  effect  kernel  function  matrices  of  equation  (4.25) 
by  one-fourth.   The  surface  over  which  the  flow  conditions 
are  specifically  satisfied  are  indicated  by  the  solid  line 
of  Figure  9b,  which  then  satisfies  the  boundary  conditions 
on  the  whole  wing/body  surface  through  symmetry.   The  body- 
wing  kernel  function  matrix  (KOTr)  follows  the  same  develop- 

BW 

ment  as  in  Section  IV. A. 3  for  the  wing  alone  and  has  the 
same  form  as  Kw  given  by  equation  (4.8). 

It  should  also  be  noted  that  for  antisymmetric  wing 
motion  the  kernel  functions  of  equations  (4.28)  and  (4.29) 
would  have  the  form 

[K]  =  [K1  -  K2  +  K3  -  K4] 


59 


^ 

/ 


«W) 


\ 


\ 


VyvaoyOcXA^leV  Flow 


\ 


^ 


y 


--"* 


A*" 

V^nq  ?o"Verrtla\ 
DouVDVe^rs 


REPKESENT/VTVON    OF   WmG/ftODY  POTENTIAL 

S\U(bA_VU\KVT\ES       , 


Special  clqAIni  Scxtvs^rvecl 


V 


9> 


t^o-Flow 
Condition  SckTv-Svyeci 


figure  qb 

WH<a/e>O0\    CslO-FLOW   BDUMO/\RY  CoMD\T\ON 


60 


Since  any  general  harmonic  motion  of  the  wing  can  be 
expressed  as  a  combination  of  symmetric  and  antisymmetric 
motion,  this  method,  as  in  the  wing  alone  case,  need  not 
be  restricted  to  the  symmetric  motion  case  considered  in 
detail  here. 

4.   Effective  Body  Length 

The  body  considered  in  this  analysis  is  idealized 
as  an  infinitely  long  cylinder  with  axis  aligned  with  the 
undisturbed  free  stream  flow.   This  is  not  an  unrealistic 
limitation,  since  in  most  practical  cases,  the  center  part 
of  a  fuselage  is  nearly  cylindrical  and  the  fineness  ratio 
of  the  fuselage  is  large  enough  so  that  the  flow  at  the 
center  part  is  nearly  the  same  as  for  an  infinitely  long 

cylindrical  body.   Experiments  in  steady  flow  have  proved 

* 

that  beyond  an  effective  body  length,  the  lift  distribution 
of  the  wing/body  combination  is  independent  of  the  body 
length  [40].   Steady  flow  analyses,  such  as  Woodward's  [4l] , 
employ  a  "wing-body  interference  region,"  where  the  body 
no-flow  boundary  condition  is  explicitly  applied  a  finite 
distance  upstream  and  downstream  from  the  wing  root  section. 

A  similar  effective  body  length  is  modeled  in  this 
approach,  where  the  body  singularity  grid  extends  from  a 
finite  distance  upstream  of  the  wing  root  to  a  point  down- 
stream of  the  effective  wake.   The  actual  body  length  which 
must  be  modeled  in  order  to  include  all  interference  effects 
is  discussed  in  Section  VI.   Since  the  kernel  function 
velocities  decrease  approximately  as  the  distance  from  the 
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singularity  cubed,  it  is  not  unreasonable  to  assume  that 
interference  effects  are  concentrated  close  to  the  wing/wake 
area. 

5.   Body  Singularities 

The  acoustic  singularity  used  to  model  the  wing  was 
a  potential  doublet,  due  to  the  requirement  for  a  pressure 
differential  across  the  wing  surface.   No  such  requirement 
exists  for  the  body  singularities,  so  that  an  harmonically 
oscillating  source  could  be  employed  as  well  as  the  dipole. 
However,  in  the  development  of  the  wing  computer  program, 
it  became  apparent  that  the  solution  to  the  matrix  equation 
(^.4)  was  very  sensitive  to  the  value  placed  on  the  kernel 
function  upwash  singularity.   Approximate  methods  used  to 
obtain  an  average  value  of  the  upwash  did  not  produce 
acceptable  loading  distributions  over  the  wing.   Only  when 
this  upwash  was  determined  numerically,  as  indicated  in 
Section  IV. A. 4  were  good  results  obtained.   Unfortunately 
a  source  singularity,  while  having  a  somewhat  simpler  form 
of  kernel  function,  can  not  be  analyzed  by  this  same  type 
of  numerical  development. 

The  body  was  therefore  modeled  by  a  network  of  harmon- 
ically oscillating  potential  doublets,  as  on  the  wing,  with 
axes  oriented  in  the  radial  direction.  Each  doublet  is  at 
the  center  of  its  control  panel,  coincident  with  a  normal- 
wash  control  point.  Determination  of  the  normalwash  (radial 
velocity  u  )  at  the  doublet  location  is  obtained,  as  with 
the  wing  singularity,  from  a  consideration  of  continuity. 
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The  total  fluid  outflow  from  the  doublet  is  u  A  ,  where 

r0  ° 
A  is  the  area  of  the  doublet  control  panel  on  the  body. 

The  amount  of  fluid  passing  back  through  the  body  at  a 
control  panel  away  from  the  doublet  is  (u  )   A   where  n 
indicates  the  axial  location  and  m  the  circumferential  loca- 
tion of  this  control  panel.   Continuity  is  then  expressed 
for  this  discrete  element  representation  as 

oo    m 

u   A   +   Z    Z   (u  )   A    =0 
r  o     ,    -    r  nm  nm 
o     n=l  m=l     nm 

where  M  equals  the  number  of  panels  located  around  the 

circumference  of  the  body.   If  the  A   ' s  are  chosen  so  that 

nm 

A   =  A   the  value  of  the  doublet  singularity  is  given  by 
nm    o  &      j     &      j 

oo    m 

u   =  -   Z    Z   (u  )  «  (4.30) 

o     n=l  m=l    r  nm 

This  continuity  equation  states,  in  effect,  that  all  the 
fluid  which  leaves  the  doublet  in  the  radial  direction  must 
pass  back  through  the  body  surface  prior  to  returning  to 
the  doublet.   As  in  the  wing  singularity  case,  the  summation 
of  normalwash  in  the  axial  direction  can  be  terminated  after 
a  reasonable  distance  in  both  the  upstream  and  downstream 
directions . 

6.   Kernel  Functions 

The  elements  of  each  of  the  kernel  function  matrices 
defined  in  Section  IV. B. 2  represent  the  normalwash  on  either 
wing  or  body  surface  due  to  a  harmonically  oscillating 
doublet  of  unit  strength  which  satisfies  the  governing 
linearized  differential  equation  (2.12).   Doublets  on  the 
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wing  are  of  course  oriented  in  the  0  (or  z)  direction, 
while  those  on  the  body  have  their  axes  in  the  radial 
direction. 

The  formulation  of  these  kernel  functions  is  developed 
in  Appendix  A.   They  are  summarized  below  with  control 
point  at  (r,9,x)  and  doublet  located  at  (r  ,9  ,x  ).   The 
body  radius  is  rR. 

(i)   Kw  -  Downwash  on  wing  due  to  unit  wing/wake  singu- 
larity. 

_r2 


Ua  = 


3  (1+  i  — ^  R)exp  {i  -£-_  [mot(x-xo)-r]  }    (4.3D 


^ttR"3       a  3  a  0 

00  oo 


where 


R 


=  -7(x"xo)2  +  e2(r2+r<;  +  2rr  )   , 


(ii)   KRW  -  Normalwash  on  body  due  to  unit  wing/wake 
singularity . 

u  =  — *-~   sin  9  (1 -r—)  (1  +  i  o  R) 

r   UttR3  n  3r      a  B2 

00 

+  (^2[rBR  i]  exP  {1  -T2  [m„(x-xo)-r]  }   (K.32) 


where 


a«,e    u  ^6 


R  -  ^(x-xq)2  +  e2(r^+r^  +  2rQrB  cos  9) 


2 

9R    p  /    —        0\ 

tt-   =  -tL-  (rn  +  r   cos  9) 
8r    r    B  T-  o 


6k 


In  both  the  above  formulations,  the  minus  sign  represents  a 
singularity  on  the  6=0  wing/wake,  and  the  plus  sign  a 
singularity  on  the  6=  tt  virtual  wing/wake. 

(iii)   KR  -  Normalwash  on  body  due  to  unit  body  singu- 
larity . 


u   --=£ 


4ttR- 


cos(e-eo)  +  |  rB  (l- cos<e-eo>)  |f 


1  +  i 


CO 


a   6' 


R 


-   <-r?)     RrB(l-cos<e-8o»  || 
a   3 


6XP    I1    -T2     fr-^-V    -    R]} 


(4.33) 


a   6 

00  H 


where 


R 


=    J(x-xo)2   +   202rg    (l-cos<9-eo» 


|f  =  ^rB(i-cos<e-eo» 


(iv)      KWR   -   Downwash   on  wing   due    to   unit    body    singu- 
larity . 
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-3 


ue  = 


4ttR- 


r   sin    eo   +   |   (rR   -    r   cos    eQ)    |§ 


1  +  i 
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V.   DESCRIPTION  OF  COMPUTER  PROGRAMS 

A.   GENERAL  DESCRIPTION 

Listings  of  the  two  computer  programs  developed  and  used 
in  this  investigation  are  reproduced  in  Appendices  D  and  E. 
The  programs  are  written  in  standard  FORTRAN  IV  language. 
Calculations  were  performed  on  the  I.B.M.  360/6 7  computer 
at  the  W.  R.  Church  Computer  Center  of  the  Naval  Postgraduate 
School.   Object  codes  were  obtained  with  the  I.B.M.  G-level 
compiler. 

Both  the  wing  and  the  wing/body  programs  are  arranged  in 
the  same  general  format.   The  MAIN  program  reads  the  input 
data  and  establishes  the  wing  or  wing/body  control  grid. 
One  or  more  subroutines  are  called  which  establish  the 
kernel  function  matrix  elements.   MAIN  calculates  the  wing 
downwash  velocity  vector  and  calls  subroutine  COMAT  to 
solve  the  matrix  downwash  equation  (4.2)  for  the  wing  velo- 
city potential  distribution.   MAIN  finally  calls  subroutine 
PRES  which  calculates  the  perturbation  pressure  distribution 
over  the  wing  by  applying  finite  difference  approximations 
to  equation  (4.5).   PRES  then  calls  subroutines  SECLM  and 
WINGLM  to  integrate  this  distribution  to  obtain  sectional 
and  wing  forces  and  moments.   These  results  are  printed  by 
PRES  and  control  is  returned  to  MAIN  for  the  reading  of 
data  cards  for  the  next  problem. 
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B.   WING  PROGRAM 

The  wing  program  computes  potential  and  pressure 
loadings  due  to  harmonic  motion  on  general  planar  wing 
configurations  from  rectangular  to  arbitrary  sweep  angles 
for  both  the  leading  and  trailing  edges.   The  program  is 
restricted,  however,  to  constant  sweep  angles,  and  to 
a  finite  tip  chord  with  minimum  taper  ratio  of  about  one- 
fifth.   Thus  delta  wing  configurations  are  excluded.   This 
restriction  is  caused  by  the  method  used  to  integrate  the 
chordwise  pressure  distribution  to  obtain  sectional  lift 
and  pitching  moment,  discussed  in  Section  IV.A.5.   The 
program  provides  for  a  maximum  of  100  control  points  on 
the  wing  and  ten  control  points  in  either  the  chordwise 

or  spanwise  directions.   Storage  requirements  are  the 

* 

equivalent  of  18,000  single  precision  complex  words,  or 
144,000  bytes  on  the  I.B.M.  360,  with  a  maximum  run  time 
of  approximately  twelve  minutes  for  one  problem.   Provision 
is  made  for  the  running  of  successive  problems  for  as  many 
sets  of  input  data  cards  as  are  provided. 

Figure  10  is  a  diagram  of  the  wing  program,  while  sub- 
routine flow  diagrams  are  presented  in  Appendix  C.   Sub- 
routine DINCO  forms  the  kernel  function  matrix  from 
equation  (4.8),  with  the  matrix  elements  defined  by  equation 
(4.3).   The  COMAT  subroutine  solves  the  linear  complex 
matrix  equation  (4.2)  by  the  Gauss-Jordan  method  with  total 
pivoting.   This  subroutine  was  written  for  systems  of  real 
equations  by  Mrs.  Sharon  Good,  David  Taylor  Model  Easin, 
as  published  in  Ref.  42,  and  modified  to  include  the  complex 
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capability  by  Mr.  Hellmut  Golde,  Department  of  Electrical 
Engineering,  University  of  Washington.   COMAT  was  further 
modified  by  the  author  for  particular  application  to  the 
wing  and  wing/body  programs. 

Input  instructions  are  presented  in  Appendix  B.   Any 
number  of  spanwise  and  chordwise  control  points,  up  to  the 
maximum  of  ten  each,  may  be  specified  for  a  semi-span  wing. 
The  wing  is  modeled  by  an  equal  number  of  chordwise  control 
points  at  each  spanwise  station.   The  program  is  limited 
by  theory  to  the  subsonic  flow  regime,  within  which  any 
mode,  amplitude,  or  frequency  of  harmonic  wing  motion  may 
be  specified,  including  the  steady  case.   Results  printed 
by  the  program  for  each  spanwise  station  are: 

(i)     Potential  doublet  strength  distribution 

4 

(ii)    Pressure  coefficient  distribution 

(iii)   Section  lift  and  pitching  moment  coefficients 

In  addition,  wing  lift  and  pitching  moment  coefficients  are 

presented. 

C.   WING/BODY  PROGRAM 

The  wing/body  program  incorporates  the  body  surface 
boundary  condition  into  the  harmonic  motion  analysis  of  a 
rectangular  wing  planform.   As  in  the  wing  program,  a  maxi- 
mum of  100  control  points  may  be  used  to  model  the  wing, 
while  up  to  130  points  are  allowed  in  the  body  control  panel 
network.   Storage  requirements  are  the  equivalent  of  ^3*500 
single  precision  complex  words,  or  3^8,000  bytes  on  the 
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I.B.M.  360.   A  maximum  run  time  of  about  20  minutes  for 
the  wing/body  problem  is  required,  while  a  wing  only  solu- 
tion (body  radius  equals  zero)  is  obtained  in  less  than 
seven  minutes.   As  in  the  wing  program,  successive  problems 
for  different  wing,  body,  and  flow  geometries  may  be  run. 

Figure  11  diagrams  the  wing/body  program,  with  indivi- 
dual subroutine  flow  charts  presented  in  Appendix  C.   The 
influence  coefficient  matrices  defined  in  Section  IV. B. 2 
are  formed  in  the  following  subroutines: 

(i)     DINCO  -  Dw  from  equation  (4.31) 

(ii)    UTHETA  -  KWR  from  equation  (4.3*0 

(iii)   URAD1  -  KBW  from  equation  (4.32) 

(iv)    URAD2  -  Kg  from  equation  (4.33) 

The  matrix  K   is  inverted  by  subroutine  COMAT  and  the  modi- 
B 

fied  kernel  function  matrix  of  equation  (4.25)  is  formed 
in  MAIN.   COMAT  is  again  called  to  solve  the  resulting  matrix 
equation  and  PRES  performs  the  same  function  with  the  same 
data  printouts  as  in  the  wing  program.   Subroutine  DINCO, 
while  performing  the  same  function  as  in  the  wing  program, 
is  here  in  a  more  simplified  form  because  of  the  restricted 
wing  geometry  of  the  wing/body  program.   This  accounts  for 
the  reduced  running  time  of  wing  only  problems  in  this 
program  as  compared  to  the  general  wing  program. 

The  wing/body  program  uses  the  same  wing  grid  and  flow 
geometry  parameters  as  the  wing  program.   Any  number  of 
body  control  points  may  be  specified  up  to  the  maximum  of 
130.   An  effective  body  lerigtn'  of  one  wing  chord  upstream 


70 


figure:  l\ 


w\vaq/^ody  program  o\/\<^a*a 


START  ^      \M9UT 
'  fc  DATA, 


± 


9>  STOP-  MO  \UPUT 


<- 


-^    UTHETA 


£~ 


^UK^OZ 


\(b 


■> 


COKA\ 


<r 


K^ 


A 

I 

N 


<r 


-^ 


D\NiCO 


<■ 


\<w 


^    eOTAATT 


<" 


<t=vc-L,w 


-e- 
ii 


\k 


p 

E 

S 


^  SECLM 


4 


->PR\NT  RESULTS 


71 


from  the  wing  leading  edge  to  one  wing  chord  downstream 
from  the  termination  of  the  effective  wake  is  modeled  by 
the  program.   Input  data  instructions  are  presented  in 
Appendix  B. 
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VI.   RESULTS  AND  DISCUSSION 

A.   WING  ANALYSIS 

1 .   Comparison  with  Lifting  Surface  Theory 

In  order  to  investigate  the  validity  of  the  discrete 
potential  element  approach  to  lifting  surface  theory,  the 
wing  program  was  run  for  a  wide  variety  of  wing/flow  geometries 
and  types  of  oscillatory  motion.   Examples  were  picked  which 
could  be  checked  against  previous  work  reported  in  the 
literature,  representing  a  variety  of  lifting  surface  theory 
approaches,  as  well  as  the  relatively  limited  amount  of 
experimental  results  which  have  been  obtained  in  the  unsteady 
field. 

4 

The  figures  presented  in  this  section  are  a  representative 
summary  of  this  investigation.   The  subsonic  flow  regime  was 
covered  from  the  incompressible  Moo=0  to  the  high  subsonic 
14^=0.9,  while  the  range  of  frequencies  varied  from  the  steady 
case  to  the  relatively  high  reduced  frequency  of  k=  1.2. 
The  wing  planforms  considered  had  a  range  of  aspect  ratios 
from  one  to  six,  sweep  angles  from  zero  to  45  degrees,  and 
taper  ratios  from  one  to  one-half.   The  types  of  symmetric 
harmonic  motion  considered  were: 

(i)    Pitching  -  rigid  body  oscillations  of  the  wing 
about  a  spanwise  axis  perpendicular  to  the  flow  direction; 

(ii)   Bending  -  oscillations  of  the  wing  as  a  beam  canti- 
levered  at  the  root  chord  in  the  first  natural  mode  of  bending. 
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(iii)   Plunging  -  rigid  body  oscillations  of  the  wing  in 
the  z  direction  at  a  mean  angle  of  attack  of  zero; 

(iv)    Flapping  -  rigid  body  roll  oscillation  of  the  wing 
simply  supported  from  a  chordwise  axis  inboard  of  the  root. 
Steady  data  were  obtained  with  the  wing  analyzed  as  a  rigid 
body  at  a  fixed  angle  of  attack.   The  results  have  been  pre- 
sented as  chordwise  pressure  distributions,  spanwise  lift 
and  pitching  moment  distributions,  and  wing  coefficients 
plotted  with  respect  to  reduced  frequency. 

In  each  of  the  figures,  the  coefficient  amplitudes  are 
normalized  with  respect  to  the  motion  as  follows: 

(i)     Pitching  -  angle  of  attack  amplitude; 

(ii)    Bending  -  wing  tip  angle  of  attack  amplitude; 

(iii)   Plunging  -  vertical  motion  amplitude; 

* 

(iv)    Flapping  -  amplitude  of  flapping  angle. 
These  represent  the  conventions  used  in  the  applicable 
references,  however  the  phase  angle  relationship  were  con- 
verted to  the  coordinate  convention  used  in  this  paper  where 
differences  occurred. 

Correlation  of  the  wing  program  results  with  existing 
lifting  surface  theory  was  in  general  quite  good.   Devia- 
tions appeared  to  come  from  the  different  methods  employed 
to  handle  the  wing  edge  singularities  in  the  pressure 
distribution.   The  lifting  surface  approaches  assume  the 
form  of  these  singularities  a  priori  in  the  choice  of  their 
loading  functions,  as  discussed  in  Section  III.B.   The 
discrete  element  approach,  on  the  other  hand,  assumes  no 
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loading  profile,  but  obtains  proper  potential  and  pressure 
distributions  through  inclusion  of  the  boundary  conditions 
in  the  problem  formulation.   As  discussed  in  Ref.  1,  much 
effort  has  been  devoted  in  lifting  surface  theory  to 
developing  the  most  numerically  efficient  and  accurate 
loading  functions,  as  well  as  to  improve  the  numerical 
methods  of  handling  the  singular  kernel  functions. 

The  first  four  figures  compare  wing  program  results  with 
one  of  the  most  recent  lifting  surface  theories.   This 
advanced  kernel  function  method,  based  on  the  acceleration 
potential,  was  presented  by  Laschka  in  19 6 3  [5]  and  further 
developed  by  Laschka  and  Schmid  for  interfering  planar 
surfaces  in  a  1967  paper  [43].   Figure  12  compares  chordwise 

pressure  distributions  on  a  swept,  tapered  wing  undergoing 

* 

pitching  oscillations  in  low  subsonic  flow  with  Laschka' s 
results  [46].   Good  correlation  is  evident  at  each  spanwise 
station,  with  the  only  variance  being  a  small  difference  in 
the  shape  of  the  pressure  distribution  near  the  leading 
edge,  as  previously  discussed. 

Figures  13  and  Ik   present  spanwise  loading  on  a  45  degree 
swept  constant  chord  wing  in  both  the  pitching  and  plunging 
modes.   Results  are  compared  with  Lashka's  work  [5]  for  both 
low  subsonic  (M  -0)  and  M  =0.8  flows.   Correlation  is 

OO  00 

again  good.   The  discrete  potential  element  approach  appears 
to  overestimate  the  amplitude  of  the  lift  and  pitching 
moment  slightly,  especially  near  the  wing  tip.   In  fact  this 
latter  trend  appears  in  almost  all  the  results  of  the  present 
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method.   The  explanation  again  would  appear  to  lie  in  the 
method  of  handling  the  wing  tip  pressure  slope  singularity, 
which  is  included  implicitly  in  the  boundary  conditions. 
Further  development  of  the  discrete  potential  element 
approach  would  definitely  have  to  include  investigation  of 
the  adequate  recognition  of  this  wing  tip  condition  in  the 
problem  formulation. 

In  Ref.  37  Laschka  compares  his  results,  for  the  wing 
configurations  of  Figures  13  and  14,  with  the  work  of  Pao 
Tan  Hsu  [4M] .   It  should  be  noted  that  the  wing  program 
results  agree  more  closely  with  Laschka' s  data,  while  fall- 
ing between  Laschka's  and  Hsu's  results  except  for  amplitudes 
at  the  wing  tip.   The  theory  here  is  further  compared  with 
experimental  results  obtained  by  Laidlaw  [45]  with  fair 
correlation. 

Figure  15  compares  wing  program  spanwise  loading  with 
Laschka's  results  from  Ref.  37  for  a  steady  swept,  tapered 
wing  in  low  subsonic  flow.   The  lift  and  pitching  moment 
coefficients  are  presented  on  an  expanded  scale,  with  the 
same  type  of  comparisons  already  noted. 

The  results  presented  in  Figures  16  through  18  are  note- 
worthy for  the  Mach  number  range  considered  (Moo=0.24  to  0.9) 
and  for  the  consistent  experimental  data  presented  from  the 
work  of  Lessing,  Troutman,  and  Menees  [47].   Spanwise  and 
chordwise  loadings  are  plotted  for  an  aspect  ratio  three 
rectangular  wing  in  the  bending  mode.   Theoretical  comparison 
is  made  to  the  lifting  surface  results  developed  by  Lessing 


76 


et  al  in  I960  from  the  original  kernel  function  methods  of 
Refs.  25  and  48. 

At  a  Mach  number  of  0.24,  the  chordwise  pressure  loadings 
obtained  from  the  two  theories  are  almost  identical.   The 
spanwise  lift  and  moment  curves  show  a  somewhat  heavier 
loading  concentration  near  the  wing  tip  for  the  discrete 
potential  element  method,  as  in  the  preceding  figures.   This 
trend  also  appears  at  Mach  numbers  of  0.7  and  0.9-   The 
M  =0.7  chordwise  pressure  distributions  show  a  difference 
in  the  representations  of  the  leading  edge  pressure  distri- 
butions outboard  of  n  -    0.5-   These  deviations  can  probably 
again  be  laid  to  the  different  methods  employed  to  handle 
the  wing  pressure  singularities.   However,  the  wing  program 
results  agree  more  closely  to  the  experimental  data  of  Ref. 
47,  than  does  the  kernel  function  approach.   No  explanation 
can  be  found  for  the  deviation  of  the  pressure  phase  angles 
in  the  wing  trailing  edge  area  for  both  the  0.24  and  0.7 
Mach  number  cases.   It  should  be  noted,  however,  that  the 
experimental  ooints  again  correlate  more  closely  with  the 
wing  program  results  in  this  area.   In  fact,  the  close  cor- 
relation of  these  experimental  data  with  the  results  of  the 
present  method  is  quite  gratifying. 

It  Is  interesting  to  observe  that  a  shock  wave  may  have 
started  to  form  on  the  wing  in  the  0.7  Mach  number  flow  down- 
stream of  the  60%   chord  line,  as  indicated  by  the  drastic  jump 
in  the  pressure  phase  angle  measurements  at  this  location  in 

Figure  17b.   The  theoretical  analyses  would  not,  of  course, 

reflect  the  existence  of  the  shock  wave. 
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The  closeness  of  the  two  theoretical  approaches  at 
M  =  0.9  is  also  noteworthy.   As  discussed  in  Section  II. C, 

oo  * 

the  linearization  of  the  basic  problem  for  perturbation 
analysis  would  appear  somewhat  questionable  this  close  to 
the  speed  of  sound.   However,  the  consistency  of  the  results 
would  seem  to  indicate  that  the  linearized  subsonic  theory 
is  still  valid  at  this  high  subsonic  Mach  number. 

Figure  19  analyzes  the  pitching  motion  of  an  aspect 
ratio  two  rectangular  wing.   Chordwise  and  spanwise  loadings 
are  compared  with  the  experimental  and  theoretical  results 
of  Laidlaw  [45].   The  latter  represents  a  numerical  treat- 
ment of  the  rectangular  wing  aspect  ratio  theory  of  Reissner 
[49,  50]  developed  in  1947.   Considering  the  fact  that 
Reissner' s  approach  basically  involves  applying  correction 

* 

factors  to  two-dimensional  results  in  order  to  account  for 
finite  aspect  ratio,  this  theory  agrees  quite  well  with  the 
discrete  potential  element  approach.   The  experimental  data, 
although  not  as  consistant  as  that  of  Lessing  [47],  agrees 
reasonably  well  with  the  Wing  Program  results.   Figure  20 
presents  spanwise  loading  for  the  same  wing  in  plunging 
motion.   Correlation  of  pressure  distributions  between 
theories  and  experimental  data  for  this  mode  of  nonsteady 
motion  is  the  same  as  for  the  pitching  case. 

The  frequency  response  of  a  rectangular  wing  in  the  flap- 
ping mode  is  presented  in  Figure  21.   Comparison  is  made  to 
experimental  and  theoretical  data  by  Woolston,  Clevenson, 
and  Leadbetter  [51] ,  who  employed  a  basic  kernel  function 
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approach.   Correlation  of  the  wing  program  results  with  the 
experimental  data  is  reasonably  good,  although  there  is  some 
difference  in  the  phase  angle  values.   However,  the  varia- 
tions of  both  lift  and  pitching  moment  with  reduced  fre- 
quency compares  very  well.   Woolston's  theoretical  points 
coincide  essentially  with  the  wing  program  curves  except 
for  the  pitching  moment  phase  angle.   No  explanation  can  be 
found  for  this  difference.   Laschka  [5]  >  in  comparing  his 
results  with  those  of  Woolston,  produced  virtually  the  same 
curves  as  the  discrete  potential  element  program. 
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2 .   Convergence 

Figures  22  through  24  present  spanwise  loading  for 
three  different  wing  planforms  as  a  function  of  the  singu- 
larity grid  density  in  order  to  show  convergence  of  the  wing 
program.   Three  types  of  motion  are  considered,  with  wing 
grids  numbering  from  25  to  100  points.   Figure  25  shows 
convergence  of  the  wing  program  for  the  steady  case  on  an 
aspect  ratio  six  rectangular  wing.   This  wing  planform  is 
used  in  the  basic  reference  for  the  wing/body  analysis  in 
steady  flow.   The  ordinate  of  these  graphs  has  been  greatly 
expanded  to  permit  a  good  evaluation  of  convergence. 

All  of  the  examples  are  seen  to  converge  at  about  the 
same  rate.   At  50  control  points  on  the  wing  the  maximum 
deviation  is  less  than  four  per  cent.   At  60  or  more  points 
the  results  are  virtually  coincident.   The  shape  of  the 
control  box  would  appear  not  to  be  a  factor,  since  varying 
the  grids  on  the  individual  wings  of  Figures  22  through  25 
also  varies  the  shape  of  the  individual  boxes.   However,  it 
would  seem  prudent  to  maintain  the  control  boxes  at  reason- 
able aspect  ratios  (less  than  about  ten)  to  assure  valid 
results.   Another  restriction  on  the  wing  grid  size  is  the 
requirement  to  integrate  the  pressure  loading  to  obtain 
section  and  wing  coefficients.   A  minimum  of  six  points  in 
either  the  chordwise  or  spanwise  directions  was  found  neces- 
sary to  adequately  define  the  loading  distribution  for  the 
integration  subroutines.   However,  the  basic  requirement  for 
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60  or  more  points  in  the  wing  grid  makes  this  latter 
restriction  a  factor  only  for  higher  aspect  ratio  wings. 

A  second  type  of  restriction  on  the  wing  grid  size  is 
caused  by  the  accuracy  of  the  inversion  method  of  subroutine 
COMAT.   Solution  of  the  downwash  matrix  equation  is  accom- 
plished by  the  Gauss-Jordan  method  employing  total  pivoting. 
With  the  relative  magnitudes  of  the  kernel  function  matrix 
elements  occurring  in  this  type  of  analysis,  the  inversion 
process  starts  to  lose  accuracy  with  a  system  of  equations 
or  order  greater  than  about  110.   Computational  accuracy 
was  verified  by  substituting  the  velocity  potential  vector 
solution  into  the  downwash  matrix  equation,  and  obtaining 
a  downwash  vector  to  compare  with  the  original  input.   Checks 
of  the  solution  for  grid  densities  up  to  150  points  were 
made,  for  various  wing/flow  conditions,  with  consistent  loss 
of  accuracy.   Thus,  the  program  has  been  limited  to  wing 
grids  of  100  control  points  or  less  to  ensure  accuracy. 
This  also  limits  the  aspect  ratio  of  the  wings  that  may  be 
analyzed;  however,  for  the  purposes  of  this  report,  aspect 
ratios  of  six  or  less  were  well  within  the  capability  of 
the  program.   If  added  capability  were  required,  a  more 
sophisticated  inversion  routine  would  have  to  be  developed. 
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3.   Wake  Effect 

Figures  26  and  27  summarize  the  effect  of  the  finite 
wake  length  considered  in  the  singularity  grid  network  as 
discussed  in  Section  IV. A. 2.   Seven  configurations  of  wing/ 
flow  geometries  were  analyzed,  covering  the  range  of  condi- 
tions considered  in  this  paper.   Wake  effects  were  included 
to  a  maximum  distance  of  ten  root  chord  lengths  downstream 
from  the  wing  trailing  edge.   In  all  cases,  wing  loading 
had  converged  to  within  one  per  cent  of  amplitude  and  one- 
half  degree  of  phase  angle  by  four  root  chord  lengths  of 
effective  wake.   This  correlates  with  Haviland's  work  [8], 
in  which  he  reported  that  the  effect  of  the  wake  on  the  wing, 
for  a  rectangular  planform  in  the  steady  case,  was  determined 
to  within  one  per  cent  in  five  chord  lengths. 

Figure  26  presents  a  specific  case  of  wing  coefficient 
variation  with  effective  wake  length.   The  maximum  deviation 
of  the  wing  coefficients,  with  respect  to  an  effective  wake 
length  of  four  root  chords,  is  shown  in  Figure  27.   From 
the  latter  figure,  it  can  be  seen  that  the  near  wake  (within 
one  root  chord  length  of  the  trailing  edge)  has  a  very  great 
effect  on  the  wing  loading.   The  wake  influence  then  dies 
out  quite  rapidly  as  the  effective  termination  point  is 
moved  downstream.   Therefore,  the  finite  wake,  taken  into 
account  in  determining  the  aerodynamic  loading  on  the  wing, 
need  only  extend  a  realtively  short  distance  downstream  in 
order  to  obtain  accurate  results. 

In  the  wing  analysis  of  this  paper,  an  effective  wake  of 
four  root  chord  lengths  was  used.   The  wing  program  will, 
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however,  incorporate  any  effective  wake  length  desired  by 
varying  an  input  parameter  (Appendix  B) .   With  a  wing  grid 
of  100  control  points,  four  root  chord  lengths  of  effective 
wake  requires  a  wake  grid  which  varies  from  400  singularities 
for  a  rectangular  wing,  to  approximately  1100  singularities 
for  a  wing  with  a  taper  ratio  of  one-fourth.   The  increased 
number  of  wake  singularities  are  required  for  a  tapered  wing 
because  of  the  decreasing  chordwise  dimension  of  the  control 
boxes  towards  the  wing  tip.   Thus  a  greater  number  of 
singularities  are  required,  at  a  spanwise  station  outboard 
of  the  root,  to  extend  four  root  chord  lengths  into  the  wake, 
than  are  required  for  a  rectangular  wing.   The  large  number 
of  wake  singularities  pose  no  computational  problem,  since 
the  effect  of  the  wake  is  incorporated  prior  to  the  solution 

4 

of  the  kernel  function  matrix  equation,  as  discussed  in 
Section  IV. A. 2. 
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B.   WING/BODY  ANALYSIS 

1.   Comparison  with  Steady  Analyses 

A  thorough  literature  search  was  made  in  an  attempt 
to  find  previous  theoretical  or  experimental  work,  in  the 
nonsteady  lifting  surface  field,  with  which  to  compare  and 
validate  the  wing/body  program  results.   Unfortunately,  no 
such  effort  which  explicitly  analyzed  body  interference 
effects  on  an  oscillating  wing's  pressure  distribution  could 
be  found.   Rodden  [38]  includes  some  body  interference  effects 
in  his  doublet  lattice  method  for  analyzing  nonplanar  con- 
figurations, by  extending  the  lifting  surface  elements  onto 
the  body  surface  near  the  wing-body  intersection.   However, 
the  actual  interference  effects  are  not  presented. 

It  was  therefore  necessary  to  compare  the  wing/body 
program  results  with  data  previously  obtained  in  the  steady 
field.   The  two  latest  works  found  were  presented  at  the 
AGARD  conference  on  Aero-dynamic  Interference  in  1970  by 
Kuchemann  [52]  and  Labrujere  [533.   Both  analyzed  an  aspect 
ratio  six  rectangular  wing  at  a  six  degree  angle  of  attack 
midmounted  on  a  cylindrical  body  of  radius  approximately 
equal  to  20%  of  the  wing  semi-span.   This  configuration  was 
also  theoretically  analyzed  by  Weber  [40]  and  experimentally 
investigated  by  Korner  [54] .   Each  of  the  theoretical  analyses 
employed  a  combination  of  surface  source  singularities  and 
vortex  line  distributions  to  construct  a  steady  kernel  func- 
tion method  solution.  - 

Figure  28  compares  the  wing/body  program  results,  for 
both  the  wing  and  wing/body  combination,  with  results 
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obtained  from  the  above  investigations.   The  discrete 
potential  element  approach  spanwise  loadings  are  seen  to 
compare  almost  exactly  with  the  theoretical  curves  obtained 
by  Weber,  while  the  experimental  data  of  Labrujere  and 
Korner  offer  good  correlation.   The  theoretical  work  of 
Labrujere  overestimates  the  spanwise  loading,  especially  in 
the  case  of  the  wing/body  combination.   The  good  agreement 
between  these  theoretical  and  experimental  results  and  the 
wing/body  program  results  provides  a  reasonable  measure  of 
confidence  in  the  discrete  potential  element  analysis  of 
this  type  of  nonplanar  configuration. 

Figure  29  indicates  the  convergence  of  the  wing/body 
program  as  the  number  of  control  points  on  the  quarter  cir- 
cumference of  the  body  was  increased  from  72  to  108.   This 
graph  is  a  greatly  expanded  representation  of  the  wing  root 
area  effects.   Outboard  of  the  25%  semi-span  point,  the 
three  grid  configurations  gave  essentially  the  same  results 
shown  in  Figure  28.   In  each  case,  body  and  wing  control 
panels  were  matched  as  closely  as  practical  in  size  and 
shape,  so  that  the  body  grid  becomes  in  essence  an  extension 
of  the  wing  network  on  the  body  surface.   Small  variations 
between  wing  and  body  panel  size  did  not  effect  the  wing 
pressure  distribution;  however,  when  wing  and  body  grids 
were  obviously  mismatched  results  became  inconsistent. 

For  all  the  data  presented  in  this  paper,  the  body  grid 
extended  from  one  chord  length  upstream  of  the  wing  leading 
edge  to  one  chord  downstream  of  the  effective  wake  termina- 
tion point.   Extending  the  body  grid  somewhat  further  in 
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either  direction  produced  no  appreciable  effect  on  the 
results.   However,  decreasing  the  length  of  the  body  grid, 
to  the  point  where  the  wing/wake  singularity  grid  extended 
beyond  the  body  grid,  caused  the  wing  pressure  distributions 
to  become  inconsistent.   This  effect  held  for  both  the 
steady  and  nonsteady  cases.   Therefore,  it  would  seem  that 
the  body  surface  which  effects  the  wing  pressure  distribu- 
tion, as  discussed  in  Section  IV. B. 3,  is  concentrated  in 
this  effective  length  between  one  wing  chord  upstream  and 
downstream  of  the  wing/wake  grid.   Further  investigation  of 
this  effect,  such  as  by  increasing  the  effective  wake  and/ 
or  the  body  grid  length  with  increased  number  of  control 
points,  was  precluded  due  to  the  accuracy  limitations  of 
the  matrix  inversion  routine  discussed  in  Section  VI. A. 2. 

4 

In  addition,  wing  body  grids  within  this  effective  body 
length  were  limited  in  the  number  of  control  points  which 
could  be  used  for  the  same  reason.   It  was,  therfore,  not 
possible  to  investigate  configurations  with  body  radii 
greater  than  about  20%  of  the  wing  semi-span.   Future  work 
with  this  method  would  definitely  require  an  improved  com- 
plex matrix  inversion  procedure  to  allow  greater  range  In 
the  wing/body  analysis. 

An  indication  of  the  effect  of  body  size,  in  the  steady 
case,  is  presented  in  Figure  30,  where  spanwise  loadings 
are  plotted,  with  a  greatly  expanded  ordinate,  for  the  wing 
previously  considered  in  combination  with  bodies  of  three 
different  radii.   It  can  be  seen  that  the  interference 
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effect  is  large  even  for  a  relatively  small  body  radius  of 
five  per  cent  of  the  wing  semi-span.   This  is  reasonable 
since  the  body  is  at  zero  angle  of  attack  and  is  providing 
no  lift  carry-through,  or  reinforcement,  for  the  wing.   This 
effect  for  bodies  of  small  radius  is  also  shown  by  Kuchemann's 
results  [52]. 
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FIGURE  28 
STEADY   WING/BOW  SPANVJ\SE  LOKD\NQ 
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2 .   Extension  to  Nonsteady  Interference 

Investigation  of  wing/body  interference  for  the  non- 
steady  case  was  conducted  for  two  modes  of  wing  motion: 
bending,  as  previously  discussed;  and  torsional  oscillations 
of  the  wing  cantilevered  at  the  root  in  an  assumed  first 
natural  mode.   In  each  case  the  coefficients  were  normalized 
with  respect  to  the  wing  tip  angle  of  attack.   It  should 
be  noted  that,  of  the  types  of  motion  normally  considered 
in  the  nonsteady  case,  these  are  the  only  two  valid  for  the 
wing/body  configuration  considered  in  this  investigation. 

Figure  31  presents  the  interference  effects  for  the  wing 
analyzed  in  Figure  16,  cantilevered  from  a  body  with  ten  per 
cent  semi-span  radius.   The  decrease  in  sectional  lift 
amplitude  follows  the  steady  case  format,  being  greatest  at 
the  root  and  essentially  disappearing  as  the  wing  tip  is 
approached.   The  change  in  phase  angle,  while  small,  remains 
relatively  constant  until  midspan  and  then  slowly  decreases. 
Figure  32  presents  spanwise  loading  in  the  root  area  for 
wing/body  configurations  with  body  radii  varying  from  five 
to  twenty  per  cent  wing  semi-span.   Again  the  trends  are 
essentially  the  same  as  in  the  steady  case,  with  a  small  but 
relatively  constant  phase  angle  difference.   Wing/body  inter- 
ference for  the  torsional  vibration  case,  with  body  radius 
equal  to  20$  of  the  wing  semi-span  is  shown  in  Figure  33- 
These  interference  effects  in  the  torsional  loading  follow 
the  same  trends  as  for  the  nonsteady  bending  mode,  as  well 
as  the  steady  case. 
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The  results  of  the  wing/body  program  show  that  nonsteady 
interference  occurs  in  the  same  way  and  via  the  same  mech- 
anism as  in  the  steady  case.   Interference  is  greatest  in 
the  root  area,  decreasing  towards  the  wing  tip.   The  dif- 
ference between  the  steady  and  nonsteady  cases  comes  from 
the  type  of  wing  loading  with  which  the  body  interfers.   For 
a  wing  at  a  steady  angle  of  attack,  the  wing  loading  is 
greatest  in  the  root  area,  and,  therefore,  the  interference 
is  of  relatively  large  magnitude.   In  the  nonsteady  case, 
the  wing  loading  is  smaller  in  the  root  area  (since  the  wing 
is  cantilevered  from  the  body  and  has  no  notion  at  the  root) 
and  increases  towards  the  tip,  as  the  wing  motion  amplitude 
increases.   Therefore,  the  interference  effects,  concentrated 
towards  the  root  area,  are  of  smaller  relative  magnitude  for 
nonsteady  motion,  as  compared  to  the  general  steady  case. 
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FVCbUKE    3\b 
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VII.   CONCLUSIONS 

In  summary,  a  discrete  potential  element  approach  to 
subsonic  numerical  lifting  surface  theory  has  been  developed 
and  shown  to  be  practical  in  predicting  the  nonsteady  load- 
ing on  harmonically  oscillating  wings.   This  approach  was 
then  extended  to  the  case  of  an  oscillating  wing,  canti- 
levered  from  a  steady  cylindrical  body,  to  investigate  inter- 
ference effects  and  show  the  versatility  of  the  basic  method. 

Correlation  of  the  wing  program  results  with  those  of 
nonsteady  lifting  surface  theory  is  generally  quite  good 
over  the  full  range  of  subsonic  flow.   Deviations,  where 
they  exist,  appear  to  come  from  the  different  methods 
employed  to  handle  the  wing  edge  singularities  in  the  pres- 
sure distribution,  which  are  assumed  a  priori  in  the  lifting 
surface  theory,  but  in  the  discrete  potential  element 
approach  are  implicit  in  the  boundary  conditions.   Primary 
deviation  appears  to  rest  in  the  handling  of  the  wing  tip 
pressure  slope  singularity,  since  sectional  lift  and  pitching 
moment  are  generally  overestimated  in  the  wing  tip  area. 
Future  work  with  the  discrete  potential  element  method  should 
include  investigation  of  a  more  adequate  acknowledgement  of 
this  boundary  condition  in  the  problem  formulation. 

A  minimum  of  about  sixty  control  points  is  required  on 
the  semi-span  wing  surface  to  achieve  convergence  of  the 
wing  program  results.   In  addition,  at  least  six  control 
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points  are  required,  in  both  the  chordwise  and  spanwise 
directions,  to  achieve  accurate  integration  of  the  pressure 
distribution  into  sectional  and  wing  coefficients.   The 
wings,  which  can  be  analyzed  efficiently  by  this  method,  are 
effectively  limited  to  medium  to  small  aspect  ratios,  because 
of  the  requirement  for  large  numbers  of  control  points,  and 
consequently  long  computer  run  times,  for  high  aspect  ratio 
wings.   The  wing  is  further  restricted  to  constant  leading 
and  trailing  edge  sweeps,  and  to  a  finite  tip  chord.   Future 
development  of  the  program  capability  should  be  pointed  at 
removing  these  restrictions. 

The  effect  of  the  wake  on  the  wing  loading  is  concen- 
trated in  the  area  just  downstream  of  the  wing.   In  all  cases 
investigated,  wing  loading  converged  to  within  one  per  cent 
in  an  effective  wake  length  of  four  root  chords  downstream 
of  the  trailing  edge,  allowing  the  wake  singularity  grid  to 
be  terminated  at  this  finite  distance  from  the  wing.   In 
addition,  the  wake  effect  is  included  in  the  wing  kernel 
function  matrix  prior  to  the  inversion  process,  through 
application  of  the  boundary  conditions  in  the  wake.   In 
this  way,  an  historical  drawback  to  the  velocity  potential 
formulation,  the  requirement  to  integrate  the  lifting 
surface  downwash  equation  over  both  wing  and  wake,  is 
removed. 

Wing/body  program  results  agree  very  well  with  both 
theoretical  and  experimental  results  for  the  steady  case. 
Results  obtained  for  nonsteady  wing  motion  appear  consistent 
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and  give  a  good  representation  of  interference  effects. 
Convergence  of  the  wing/body  program  is  achieved  with  a  body 
control  panel  grid  of  essentially  the  same  format  as  the 
wing  grid. 

The  effective  body  length  which  causes  interference  in 
the  wing  pressure  loading  extends  approximately  from  one 
chord  length  upstream  of  the  wing  leading  edge  to  one  chord 
length  downstream  of  the  effective  wake  termination  point. 
Increasing  the  body  length  modeled  did  not  effect  the 
pressure  distribution,  while  decreasing  this  length  so 
that  the  body  grid  did  not  extend  beyond  the  wing/wake 
singularity  grid  in  either  the  upstream  or  downstream  direc- 
tions caused  inconsistent  results  to  occur.   This  analysis 
was  somewhat  limited  due  to  numerical  restrictions  on  the 

« 

allowable  size  of  the  body  grid. 

In  the  numerical  analysis,  the  number  of  control  points 
on  the  wing  and  the  body  was  limited  due  to  loss  of  compu- 
tational accuracy  in  the  matrix  inversion  routine  of  sub- 
routine COMAT.   This  restricted  the  scope  of  the  wing/body 
interference  investigation  both  as  to  the  length  and  radius 
of  body  which  could  be  considered.   Accurate  solution  of 
the  wing  downwash  equation  and  inversion  of  the  body  influence 
coefficient  matrix,  is  limited  to  coefficient  matrices  of 
order  less  than  110.   Future  wing/body  analyses  by  this 
method  would  definitely  require  a  more  sophisticated  inver- 
sion routine  capable  of  handling  much  larger  body  grids. 
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Within  the  limitations  noted,  the  wing/body  program 
results  show  that  interference  effects  are  significant,  and 
follow  the  same  format,  in  both  the  steady  and  nonsteady 
cases.   The  decrease  in  wing  pressure  loading  amplitude, 
caused  by  the  body's  presence,  is  greatest  in  the  root  area, 
decreasing  towards  the  tip.   Differences  in  phase  angle, 
while  small,  exist  over  the  entire  wing.   Therefore,  these 
interference  effects  are  of  importance  to  three-dimensional 
analyses  of  wing/body  configurations. 

Future  development  of  the  discrete  potential  element 
approach  should  include  the  analysis  of  oscillating  wings 
with  control  surfaces.   This  approach  appears  to  provide  a 
direct  means  of  including  the  control  surface  boundary  con- 
ditions, without  the  problems  associated  with  properly 

4 

defining  singularity  effects  in  the  continuous  loading 
formulation. 
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APPENDIX  A 
KERNEL  FUNCTION  DEVELOPMENT 

For  the  purposes  of  nonsteady  lifting  surface  analysis, 
a  kernel  function  is  defined  as  the  normalwash  on  a  surface 
due  to  a  unit  oscillatory  acoustic  singularity  which 
satisfies  the  linearized  convective  wave  equation  (2.9). 
The  singularities  may  be  elementry  radiators  of  zero  order 
representing  simple  point  sources,  commonly  used  in  super- 
sonic analyses,  or  first  order  radiators,  namely  dipoles  or 
doublets,  whose  axes  are  normal  to  the  surface,  used  in 
subsonic  analyses  as  in  this  paper.   Development  of  the 
general  form  of  the  doublet  kernel  function  was  first 
accomplished  by  Kussner  in  19^0  [2].   The  development  in 
this  appendix  follows  Kussner 's  method  and  specializes  the 
results  to  the  forms  of  the  kernel  function  employed  in 
the  discrete  potential  element  approach  for  both  the  wing 
and  body  surfaces. 

In  summary,  the  solution  to  the  linearized  wave  equation 
for  a  stationary  acoustic  singularity  is  first  extended  to 
the  case  of  a  moving  singularity  in  a  space  fixed  coordinate 
system  through  application  of  a  Lorentz  coordinate  transfor- 
mation invariant  with  respect  to  the  speed  of  sound.   This 
solution  is  then  transferred  to  the  moving  wing  (or  body) 
fixed  coordinate  system  through  a  Galilean  transformation. 
The  results  are  subsequently  reduced  for  the  harmonic  case 
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to  the  forms  of  the  kernel  functions  necessary  in  the  analy- 
sis discussed  in  Section  IV. 

1.   Stationary  Singularity 

The  perturbation  analysis  of  Section  II  led  to  the 
linearized  convective  wave  equation  (2.9)  ,    which  for  the 
stationary  case  (U   =0)  has  the  familiar  form 
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Since  the  pressure  and  velocity  potential  are  linearly 
related  by 
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the  purturbation  pressure  also  satisfies  the  wave  equation 
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Considering  a  stationary  source  at  the  origin  of  the  coor- 
dinate system,  equation  (A. 3)  becomes 
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where 
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This  equation  has  the  well  known  solution 
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where  F  and  G  are  arbitrary  functions.  Prom  the  boundary 
conditions  of  the  problem,  only  the  solution  representing 
outgoing  waves  (F)  is  admitted. 

Representing  the  oscillating  source  as  a  sphere  of  radius 
r,  expanding  and  contracting  with  a  rate  of  flow  away  from 
the  surface  defined  by 

S(t)  =  ^Trr^  u  (t) 
b   r 

The  momentum  equation  (2.2)  becomes  in  linearized  form 

r   3r  2  P°°      3t 

r 


or 

i|z  .4 !!=«     r  (A-6) 

r    9r  2  ji,2    dt  b 

r  ^r, 

b 

If   r,     is    very    small,    — p    is   much    larger   than  -^  at    r=rfe s    and 

equation    (A. 6)    can   be   taken   as 
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from  equation  (A. 2).   Equations  (A. 7)  represent  the  pressure 
and  velocity  potential  distributions  of  a  stationary  oscil- 
latory source  of  strength  S. 

2.   Moving  Singularity 

Consider  a  source  moving  with  uniform  velocity  U  ,  with 
respect  to  the  surrounding  fluid,  in  the  direction  of  the 
negative  x  axis.   If  Q(x,y,z,t)  is  the  source  distribution 
density,  the  continuity  equation  (2.1)  becomes 

||+  VpU  =   Q 

and  the  linearized  wave  equation  (A. 3)  becomes 

2- 
n2-    1  9  p      3Q  /.  o\ 

v  p  "  —  772  =  "  at  (A'8) 

a   dt 

00 
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This  source  distribution  can  be  represented  by 
Q(x,y,z,t)  =  pJS(t)  6(x+Uoot)  6(y)  6(z) 

where  6  represents  the  familiar  Dirac  delta  function.   Making 
use  of  the  linear  relationship  between  p  and  <j>  in  equation 
(A. 2),  equation  (A. 8)  may  then  be  expressed  in  terms  of  the 
velocity  potential  as 

2- 

V2$  -  K  ^-4   =  S(t)  6(x+Uoot)  5(y)  6(z)         (A. 9) 
a  dtd 

oo 

A  Lorentz  transformation,  scaled  with  respect  to  the  compres- 
sibility factor  S,  can  now  be  used  to  reduce  the  right  hand 
side  of  equation  (A. 9)  to  that  of  a  stationary  source.   The 
transformed  coordinate  system  (primed)  is  defined  as  follows: 
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x'  =  -\    (x+U  t) 
B^ 

yT  =  y/B 

(A. 10) 
z'  =  z/B 

i      u 
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Performing  this  transformation,  equation  (A. 9)  takes  the 

form 

2- 
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Equation  (A. 9)  has  therefore  been  reduced  to  the  form  of  the 
wave  equation  for  a  stationary  source  in  the  primed  coordi- 
nate system,  which  from  equation  (A. 7)  has  the  solution 

$(r',t')  =  ^?-  S(t'  -  r'/aj  (A. 12) 

In   terms    of  the    space    fixed   coordinate    system,    the    variables 
of  equation    (A. 12)    are 

r'   =  ^2  V(x+urot)2  +   B2(y2+z2) 

-i  u 

t'    =   ~    (t    +   -p    x) 

B  a 

h  00 

Transferring  the  solution  (A. 12)  to  the  wing  fixed  coordinate 

system  now  requires  the  further  Galilean  transformation 

x  =  x+U  t,  which  produces  the  final  result 

*(r,t)  =  t^j  s[t+-i-2  (M„x-r)]  (A.  13) 

a.  B 


13^ 


where 

r   =  ~\Jx2   +     62(y2  +  z2) 

All  variables  are  referenced  to  the  wing  fixed  coordinate 
system. 

3.   Wing  Fixed  Periodic  Singularities 

If  an  harmonically  oscillating  source  of  unit  strength 
(S  =  le    )  is  located  at  the  origin  of  the  wing  fixed 
coordinate  system,  the  velocity  potential  distribution  given 
by  equation  (A. 13)  is 

?s(x,y,z3t)  =  -  -^  exp  {ico  [t  +  -^  (M^x-r  )]}  (A.  lH) 

Considering  now  only  the  spatial  variation  of  the  velocity 

potential  of  a  unit  source  located  at  a  general  point 

(x  .y  ,z  )  (time  dependence  e  '   having  been  factored  out) 
000 

4>s(x,yaz)    =   -j=±  exp    {l   -^2    [mJx-x0>  "  R]}  (A-15) 

a  p 
00 

where 

R  =  ^/(x-x^2  +  62[(y-y0)2  +  (z-zQ)2] 

The  potential  distribution  of  a  dipole  or  doublet  is 
related  to  that  of  a  source  by 

*d  =  "8rT 

where  n  is  the  direction  of  the  axis  of  the  dipole.  There- 
fore from  equation  (A. 15),  the  velocity  potential  of  a  unit 
dipole  oriented  in  the  z  direction  is 
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a*s         63 

*d(x'y'z)=-§F-=r73  (z-zo} 

O         47Tk 


to 


1  +   i- 

\  a   3 


aRjexpfi^f 


Moo(x-xq)  -  I 


(A. 16) 


The  downwash  at  a  point  on  the  wing  (x,y,o)  from  such  a 
doublet  located  in  the  plane  of  the  wing  (x  ,y  ,o)  is  given 


by 


w(x,y,o)    = 


9<J> 


9z 


z=0,    z    =0 

'       o 


or 


w  =  ^      1+  i-^R 
^ttR^    I  a    3 

oo 


(i-^T    [Moo(^-^0)-P]}       (A. 17) 


expii ^ 


a   3 

00  M 


where 


R  =  -J(x-x)2  +   32(y-y    )2 


This  is  the  standard  form  of  the  kernel  function  in  cartesian 
coordinates  normally  used  in  the  velocity  potential  formula- 
tion of  lifting  surface  theory. 

k.      Kernel  Function  for  Wing/Body  Doublets 

To  analyze  the  wing/body  problem,  a  transformation  from 
cartesian  to  cylindrical  coordinates  was  made  as  discussed 
in  Section  IV. B  and  shown  in  Figure  8.   Thus 

x  =  x 

y  =  r  cos  0 
z  =  r  sin  0 

relate  the  equations  of  the  previous  section  to  the  (r,0,x) 
system.   A  unit  doublet  with  axis  in  the  z  direction  located 
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at  (r0>90>x0)  nas  a  potential  distribution,  from  equation 
(A. 16),  in  the  new  coordinate  system  of 


<j>d(r,6,x)  =  -£— =■  (r  sin  6-  rQsin  0q)  1  +  i  — ^ R 


4ttR- 


a  B' 


where 


exp  1  1 


I1  -^2  [M-<x-xo>-R]}         (A'18) 


a„6 


R  =  -\/(x-x  )2  +  32(r2+r2  -  2rr   cos  <6-6  >) 


For  a  doublet  located  on  the  wing  at  (r  ,  6  =  0  or  tt,x  ),  the 
downwash  on  the  wing  at  (r,o,x)  is  given  by 


urt  =  — 


8*d 


i  =  0,  0  =0  or  tt 
o 


0    r  90 
which  applied  to  equation  (A. 18)  results  in 


u 


_    -6' 


6         4,R3 


1  + 
\       a  e 


i^2    RJexpji   -^    [Mro(x-xo)-  Ft]}  (A. 19) 


where 


R  =  sj(x-x    )2    +  62(r2+r2  +  2rrQ) 


-  for  0=0 
+  for  0  =  TT 

This  is  of  course  the  same  as  equation  (A. 17)  in  cartesian 
coordinates  except  for  R  which  indicates  the  coordinate 
transformation. 
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The  normalwash  on  a  cylindrical  body  of  radius  rv, , 
coaxial  with  the  x  axis,  at  a  control  point,  (rR,9,x),  due 
to  the  wing  doublet  is  given  by 


a+; 


u      = 


r        dr 


r=rD,    6    =0    or   tt 
B'      o 


or 


where 


-6" 
u      =   _ 

r  4TTR3 


in    6      1- 


3r 


Wl+i   -^R 


R      8r 


*„* 


to 


r^R 


a„e 


2      XB 


|a]exp{i  -«     [m(x-x0)-r]J 

-J  a   p 

00 


(A. 20) 


R   =    V(x-x    )2   +    B2(r2+r2   +    2rDr„    cos    0) 


Bo 


Bo 


tt—  =  -ft    (rD    +    n      cos    6) 
9r       R         B  o 


-    for    9=0 

+    for    0    =    it 


To  determine  the  body  singularity  kernel  function  forms, 
the  velocity  potential  of  a  source  is  expressed  in  cylindri- 
cal coordinates.   From  equation  (A. 15) 

*s(r>G>x)  =  M  exp  I1  ~^2    L^Jx~xo)    ~   R]f    (A'21) 

a  3 


where 


R 


=  -V(x-x  )2  +  62(r2+r2-2rr   cos  <0-0  >) 
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The  potential  distribution  of  a  doublet  oriented  in  the 


radial  direction  and  located  at  (r  ,9  ,x  )  is  obtained  from 

o      o      o 


8<t> 


<J>.    = 


d        8r 


which   applied   to   equation    (A. 21)    gives 


4>,(r,9,x)    =   -fi-s-    [r  -r   cos    (6-8 jl    1+1-~rJ 


exp  {i   —^-2    [m(x-xq)  -  rJ} 


(A. 22) 


a„8 


with  R  as  above.   From  this  potential  formulation,  the 
following  kernel  functions  can  be  obtained  for  a  unit  doublet 
with  axis  in  the  radial  direction,  located  on  the  body  sur- 
face at  (rB,  6q,  xq). 

The  downwash  on  the  wing  at  (r,o,x)  is  given  by 


u6         r    39 


i  =  0,    r   =rn 
5       o      B 


or 


where 


u, 


r 


4ttR- 


r   sin    9     +  ^    (vr    cos 
On  d 


,  )iRl 

o;39j 


[ 


1+  i 


CO 


R 


a   B' 


ui 


2 


exp  \  l 


w 


a   6' 


a  8 

OO  / 

[m(x-xq)  -  r] 


R    (rB-r    cos    9q)    |f 


(A. 23) 


R   =  -^(x-y :    )d    +    82(r2  +  rg-2rrB    cos    9q) 
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8R   _   -B 

ae   -  —  rrB  sin   9o 


The   normalwash   on   the   body    at    (rR,6,x)    is    given   by 
8<t> 


or 


where 


B 


u     = 


r        8r 


r=rB>    ro=rB 


u 


h  77  »]  "  I 


cos(e-e  )  +  |  rR(i-cos  <e-e  >)  — 

o         R      B  o         3r 


] 


a   6 


*M   RrR(l-cos    <e-e    >)   |£ 
2  /        B  oar 


exp    j  i 


to 


a   6' 

00 


[Moo(x-xo)    -    r] 


(A. 24) 


R      =   ~J(x-x    )2    +    232r^    (1-cos    <6-6^>) 


B 


o 


|£  =  £-  rn(i-cos  <e-e  >) 

9r         r       B  o 


Equations  (A. 19),  (A. 20),  (A. 23),  (A. 24)  represent  the 
kernel  function  formulations  used  in  the  discrete  potential 
element  analysis  of  the  wing  and  wing/body  programs. 
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APPENDIX  B 
INPUT  INSTRUCTIONS  FOR  COMPUTER  PROGRAMS 

Instructions  for  preparing  input  data  for  both  the  wing 
and  wing/body  programs  are  presented  in  this  appendix.   The 
field  location  and  format  for  each  input  quantity  is  speci- 
fied.  Any  set  of  units  may  be  used  for  geometric  dimensions, 
displacements,  and  acoustic  velocity  as  long  as  they  are 
consistent.   Any  number  of  problems  may  be  run  in  sequence, 
with  the  programs  terminating  when  a  new  data  set  is  not 
available  to  be  read. 

1.   Wing  Program 

a.   First  card:   FORMAT (3P10 . 4 ,2110 ,F10\ 4 ) 

Column   1-10   11-20   21-30    40    50   51-60 
Name      XM     F      A     ICH   JCH    ALPH 

XM      Mach  number. 

F       Circular  frequency  (rad/sec). 

A       Speed  of  Sound. 

ICH     Indicator  for  second  card  input  data. 

=  0  for  new  data;  program  reads  second  card. 
=  1  for  same  data  as  previous  problem;  no 
second  card  required. 

JCH     Indicator  for  wing  displacement  input  data. 
=  0  for  new  data;  program  reads  third  and 
subsequent  wing  displacement  cards. 
=  1  for  same  data  as  previous  problem  or 
for  pitching  motion;  no  wing  displacement 
cards  required. 

ALPH    Amplitude  of  pure  pitching  motion. 

>  0  for  pitching  motion;  program  computes 
downwash  for  pitching  mode  (JCH  must  equal  1) 
=  0.0  for  general  wing  motion;  program  com- 
putes downwash  from  wing  displacement  input. 
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b.   Second  card:   FORMAT( 415 ,5F10 . 4 ) 

Column   1-5  6-10  11-15  16-20  21-30  31-40  41-50  51-60  61-70 
Name      N    M    MB    NC     XC    XB     BR    Gil    G22 

N       Number  of  chordwlse  control  points  (maximum  ten). 

M       Number  of  spanwise  control  points  (maximum  ten). 

MB      =0  for  wing  program. 

NC      Significant  wake  length  downstream  from  wing 
trailing  edge  (root  chord  lengths). 

XC      Wing  root  chord. 

BR      =0  for  wing  program 

XB      Wing  semi-span. 

Gil     Leading  edge  sweep  angle  (deg.). 

G22     Trailing  edge  sweep  angle  (deg.). 

Gil  and  G22  positive  for  downstream  sweep. 


c.   Third  and  subsequent  wing  displacement  cards: 
FORMAT (7F10. 4) 

Column   1-10  11-20  21-30  31-40  41-50  51-60  61-70 

Name    Z(1)Z(2)  Z(3)  Z(4)  Z(5)   Z(6)   Z(7) 

Z(8)  Z(9)  Z(10)  fourth  card 

Z(ll)  Z(12)  Z(13)  etc.  fifth   card 

Displacements  (wing  motion  amplitudes)  are  read  span- 
wise  starting  from  the  wing  root  at  the  leading  edge 
A  new  card  must  be  started  for  each  spanline  of  data 
(example  above  is  for  ten  spanwise  points). 


2.   Wing/Body  Program 

a.  First  card:   same  as  wing  program. 

b.  Second  card:   FORMAT ( 415 , 3F10 ,15 ) 

Column   1-5  6-10  11-15  16-20  21-30  31-40  41-50  51-55 
Name     N    M    MB    NC     XC     XB     BR    MS 

N       Number  of  wing  chordwise  control  points 
(maximum  ten) . 

M       Number  of  wing  spanwise  control  points 
(maximum  ten) . 
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MB      Total  number  of  control  points  on  body 

quarter-circumference  surface  (maximum  130) 

NC      Significant  wake  length  downstream  from  wing 
trailing  edge  (chord  lengths). 

XC      Wing  chord. 

XB      Wing  semi-span. 

BR      Body  radius. 

MS      Number  of  control  points  along  body  quarter 
circumference . 


c.   Third  and  subsequent  wing  displacement  cards:  same 
as  wing  program. 
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APPENDIX  C 
COMPUTER  PROGRAM  FLOW  DIAGRAMS 
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APPENDIX  D 
LISTING  OF  THE  FORTRAN  CODE  FOR  THE  WING  PROGRAM 
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ABSTRACT 
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monically oscillating,  medium  to  low  aspect  ratio  wings, 
unique  method  of  including  the  wake  effect  in  the  wing 
kernel  function  matrix  prior  to  solution  of  the  singular 
integral  downwash  equation  was  devised,  thus  greatly 
simplifying  the  velocity  potential  formulation.   In  addi- 
tion, termination  of  the  effective  wake  a  finite  distance 
downstream  of  the  wing  was  investigated,  with  wing  loading 
found  to  converge  to  within  one  per  cent  in  an  effective 
wake  length  of  four  root  chords. 

This  discrete  element  method  has  also  been  extended  to 
the  case  of  an  oscillating  wing,  cantilevered  from  a  cylin- 
drical fuselage,  to  investigate  nonplanar  interference 
effects.   This  interference  in  wing  loading,  while  of  rela- 
tively small  magnitude,  does  exist  in  both  pressure  amplitude 
and  phase  angle  distributions,  and  is,  therefore,  of  impor- 
tance in  three-dimensional  stability  analyses  of  wing/body 
configurations . 
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